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1. 


INTRODUCTION 


The  objective  of  the  project  is  to  provide  a  better 
understanding  of  the  interactions  among  the  chemical,  radi¬ 
ative  and  dynamical  processes  which  govern  the  behavior  of 
trace  gases  in  the  stratosphere  and  troposphere.  The  work 
involved  both  model  development  to  incorporate  the  proper 
mechanisms  as  well  as  the  interpretation  of  model  results 
through  the  comparison  with  available  observations. 

The  report  is  organized  in  seven  sections  with  two 
appendices.  Section  2  provides  a  general  review  of  numer¬ 
ical  modeling  of  tropospheric  and  stratospheric  processes 
with  emphasis  on  the  more  recent  theories  and  experimental 
data.  Section  3  gives  a  detailed  description  of  the  AER 
one-dimensional  photochemical  diffusion  model  which  includes 
the  fully  time-dependent  diurnal  code  developed  under  the 
present  contract.  Section  4  compares  the  calculated  one¬ 
dimensional  results  with  available  observations.  Section  5 
presents  a  comprehensive  study  of  the  tropospheric  and 
stratospheric  sulfur  budget  with  emphasis  on  the  chemical 
interaction  of  reduced  sulfur  compounds  with  Oil  radicals. 
Section  6  focuses  on  the  possible  ozone  perturbation  due  to 
man's  activities  including  atmospheric  release  of  fluoro¬ 
carbons,  brominated  halocarbons  and  nitrogen  oxides.  Finally, 
in  Section  7,  the  basic  approach  for  two-dimensional  zonal- 
mean  models  is  reviewed.  In  addition,  the  development  of 
the  AER  2-D  model  is  presented  together  with  the  preliminary 


model  results  simulating  the  present  atmosphere. 

Appendix  A  and  Appendix  B  contain  descriptions  of  tin 
AER  1-D  and  2-D  models  respectively .  The  results  from  tin: 
2-D  modeling  efforts  arc  being  prepared  for  submission  to 
technical  journals  for  publication. 


2. 


ATMOSPHERIC  MODELING: 


GENERAL  APPROACH 


2. 1  Background 

Atmospheric  models  are  designed  to  simulate  the  behav¬ 
ior  of  atmospheric  gases  by  solving  the  system  of  mass  con¬ 
tinuity  equations 


9n. 

i 

at 


-  V  •  <t> .  +  Q  . 
~  ~l  1 


(2.1-1) 


for  i  =  1,  .  .  .  ,  N  where  N  is  the  total  number  of  species 
considered,  n^  is  the  number  density  of  the  i  species,  'K 
is  the  flux  of  the  species  due  to  atmospheric  motion,  is 
the  net  chemical  production  or  loss  rate.  The  system  of 
equations  is  coupled  because  the  term  for  species  i  will 
depend  on  other  species  concentrations.  Equation  (2.1-1)  can 
be  interpreted  as  giving  the  time  rate  of  change  of  the  local 
species  concentration  in  terms  of  two  components,  the  part 
due  to  chemical  interactions  and  V  •  the  part  due  to  atmo¬ 
spheric  motions. 

Independent  of  the  dimensionality  of  the  model  and  the 
numerical  method  used,  a  chemical  scheme  must  be  set  up  to 
define  the  term  Q^.  The  chemical  scheme  will  specify  the 
number  of  species  considered  in  the  computation  and  also 
lists  the  reactions  among  the  various  species.  In  addition, 
boundary  values  have  to  be  specified  for  the  differential 
equations.  The  setting  up  of  the  scheme  is  guided  by  atmo¬ 
spheric  gas  concentration  measurements  as  well  as  kinetic 


measurements  in  the  laboratory.  In  Section  2.2  the  atmo¬ 
spheric  chemistry  will  be  reviewed  with  emphasis  on  recent 
developments . 

The  description  of  the  flux  term  depends  on  the  dimen¬ 
sionality  of  the  model.  In  Section  2.3,  we  will  briefly 
discuss  the  various  possible  approaches.  The  exact  approach 
and  the  numerical  methods  currently  employed  in  the  AER  1-D 
and  2-D  models  will  be  discussed  in  Sections  3  and  7  respec¬ 
tively  . 


2 . 2  Review  of  Atmospheric  Chemistry 

Stratospheric  chemistry  has  been  reviewed  extensively 
by  several  authors  (Nicolet,  1975;  Johnston  and  Pokolske, 
1978;  Logan  et  al . ,  1978;  Thrush,  1978).  An  excellent  sum¬ 
mary  of  information  available  up  to  the  mid-1979  is  given 
by  NASA  (1979) .  We  shall  focus  here  on  subsequent  develop¬ 
ment,  with  particular  emphasis  on  the  coupling  aspects  of 
various  chemical  processes  which  may  affect  the  stratospheric 
ozone  distributions. 

The  first  theoretical  model  for  stratospheric  ozone 
was  put  forward  by  Chapman  (1930)  who  invoked  four  reactions 
involving  pure  oxygen  chemistry: 


O  2  +  h  -*•  0  +  0 

o  +  o2  +  m  -►  o3  +  m 

03  +  hv  +  02  +  0 
0  +  03  ►  20 2 


(2.2-1) 

(2.2-2) 

(2.2-3) 

(2.2-4) 


Assuming  photochemical  equilibrium,  the  concentration  of 
based  on  reactions  (2.2-1)  -  (2.2-4)  is  given  by 

[03]  =  [02]2  [M]}*5  '  (2.2-5) 

where  and  are  the  photolysis  rates  for  reactions  (2.2-1) 
and  (2.2-3)  respectively,  and  k2  and  the  rate  constants 
for  reactions  (2.2-2)  and  (2.2-4)  respectively. 

Chapman's  simple  model  successfully  predicts  the  approx¬ 
imate  altitude  (^25  km)  at  which  the  ozone  peaks  and  gives 
an  ozone  overburden  about  a  factor  of  2  larger  than  observa¬ 
tions.  The  major  deficiencies  of  Chapmans'  model  are  that: 

i)  it  fails  to  account  for  the  observed  latitudinal 
and  seasonal  variations,  and 

ii)  it  tends  to  underestimate  the  sinks  for  O^. 

The  first  deficiency  is  thought  to  be  caused  by  the  neglect 
of  dynamical  transport  in  Chapman's  model,  while  the  second 
is  due  to  the  omissions  of  important  chemical  processes 
which  involve  hydrogen,  nitrogen,  chlorine,  and  perhaps 
bromine  species. 

The  importance  of  hydrogen  radicals  (H,  OH,  110 ^ )  in 
ozone  chemistry  was  recognized  by  Bates  and  Nicolet  (1950). 

It  is  now  believed  that  odd  oxygen  (0,  O^)  may  be  removed 
by  reactions  such  as 


H  +  03  -*■  OH  +  02  (2.2-6) 

OH  +  0  H  +  02  (2.2-7) 

2-3 


(2.2-8) 


H  +  02  +  M  ->•  H02  +  M 

H02  +  0  ->  OH  +  02  (2.2-9) 

OH  +  03  ->  H02  +  02  (2.2-10) 

H02  +  03  -*  OH  +  202  .  (2.2-11) 

Reactions  (2.2-10)  -  (2.2-11)  are  important  only  in  the 
lower  stratosphere  (z  <  28  km)  while  reactions  (2.2-6)  - 
(2.2-9)  are  the  major  sinks  for  ozone  above  40  km. 

The  source  for  atmospheric  hydrogen  radicals  (HOx)  is 
mainly  due  to  the  reaction  (Cadle,  1964;  Hampson,  1964) 

0(1D)  +  H20  ->  20H  (2.2-12) 

where  O^D)  is  formed  by  photolysis  of  03  at  wavelengths 

O 

shorter  than  3100  A.  The  atmospheric  sinks  for  HO^  include 
the  following  reactions: 


(2.2-13) 
(2.2-14) 
(2.2-15) 
(2.2-16) 

while  reaction  (2,2-13)  is  important  throughout  the  strato¬ 
sphere,  the  contributions  of  reactions  (2.2-14)  to  (2.2-16) 
as  sinks  for  HOx  are  restricted  mainly  to  the  lower  strato¬ 
sphere  (15-30  km) . 


OH  +  H02  •*  H20  +  02 

OH  +  HN03  -*  H20  +  N03 

OH  +  H202  >  H20  +  H02 

oh  +  ho2no2  -+  h2o  +  no2  +  o2 


Recent  laboratory  work  seems  to  indicate  that  rate 
constants  for  reactions  (2.2-13)  to  (2.2-15)  may  be  signif¬ 
icantly  faster  than  the  previous  accepted  values  (cf.  NASA, 
1979).  We  summarize  in  Table  2.2-1  the  more  recent  rate 
data  for  reactions  (2.2-13)  to  (2.2-15)  along  with  NASA's 
(1979)  recommendations. 

The  rate  constant  k^g  for  reaction  (2.2-16)  has  not 
been  measured,  although  a  tentative  upper  limit  of  about 
3x10  ^  cm3  s  1  was  set  by  Trevor  et  al.  (1980)  .  Recent 
calculations  by  Sze  and  Ko  (1980c)  show  that  reaction  (2.2-16) 

could  be  an  important  sink  for  H0x  if  k^g  exceeds  about 

.  in-12  3  -1 

1x10  cm  s 

Other  major  catalytic  cycles  for  removal  of  atmospheric 
ozone  include  the  sequences 


NO  +  03 

->• 

N02  +  02 

(2.2-17) 

no2  +  0 

-> 

NO  +  02 

(2.2-18) 

Cl  +  o3 

-b 

CIO  +  o2 

(2.2-19) 

CIO  +  0 

-> 

Cl  +  o2  . 

(2.2-20) 

The  source  for  stratospheric  N0x  (NO,  N02,  HN03,  N03, 
N2°5^  mainly  ^ue  to  t^ie  reaction  (Nicolet  and  Vergison, 
1971;  McElroy  and  McConnell,  1971;  Crutzen,  1971) , 


Thu  precursors  to 


with  N20  emanating  from  the  biosphere, 
stratospheric  chlorine  are  a  variety  of  halocarbons  (CH^Cl, 

CC14,  CH3CC13.  CFC13,  CF2C12)  as  suggested  by  Rowland  and 
Molina  (1975) . 

The  cycles  of  hydrogen,  nitrogen  and  chlorine  species 
are  strongly  coupled.  For  instance,  the  reaction 

OH  +  N02  +  M  -*■  HN03  +  M  (2.2-22) 

ties  up  active  N02  to  form  the  relatively  inert  species 
HN03  which  does  not  directly  attack  ozone.  On  the  other 
hand,  the  reaction 

OH  +  HC1  -*•  H20  +  Cl  (2.2-23) 

releases  active  chlorine  (Cl)  from  the  inert  form  HC1 .  Thus, 
an  increase  in  OH  tends  to  reduce  the  N0X  cycle  [reactions 
(2.2-17)  +  (2.2-18)]  but  augments  the  C1X  cycle  [reactions 
(2.2-19)  +  (2.2-20)]  as  sinks  for  ozone. 

The  coupling  of  NO  with  H02  and  CIO  through  the  reactions 

NO  +  H02  ■+  N02  +  OH  (2.2-24) 

and 

NO  +  CIO  +  Cl  +  N02  (2.2-25) 

is  of  considerable  importance  in  the  study  of  03  perturbations 
associated  with  increases  in  NOx  and  C1X.  Note  that  the  fol¬ 
lowing  sequence 
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(2.2-10) 


OH  +  03 

> 

no2 

+  O 

H02  +  NO 

-► 

OH  + 

NO 

N02  +  hv 

'► 

NO  + 

0 

+  °2  +  M 

-► 

°3  + 

M 

does  not  lead  to  removal  of  03. 
conserved  through  the  sequence 

Cl  +  03  CIO  +  o2 
CIO  +  NO  -  Cl  +  N02 


(2.2-24) 

(2.2-26) 

(2.2-3) 

Similarly,  odd  oxygen  is 

(2.2-19) 

(2.2-25) 


followed  by  (2.2-26)  and  (2.2-3).  Thus,  an  increase  in  NO 
may  either  reduce  O^  through  reactions  (2.2-17)  and  (2.2-18) 
or  increase  ozone  through  the  coupling  of  NO  with  H02  (reac¬ 
tion  (2.2-24)]  and  CIO  [reaction  (2.2-25)]. 

The  reactions  of  CIO  with  H02  and  N02, 


and 


CIO  +  H02  -*•  H0C1  +  02 

cio  +  no2  +  m  >  ciono2  +  M 


(2.2-27) 

(2.2-28) 


introduce  two  other  examples  of  coupling  between  ClX,  HOx  and 

N°x  cycles.  Note  that  reaction  (2.2-27)  followed  by  the 
sequence 


H0C1  +  hv 

ci  +  o3 

OH  +  03 


-*•  OH  +  Cl 
♦  CIO  t  o2 
'*  H02  +  02 


(2.2-29) 

(2.2-19) 

(2.2-10) 
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is  equivalent,  to 


20  j  *  30  2  • 

The  importance  of  reaction  (2.2-28)  involves  the  removal 
of  two  radicals  CIO  and  N02  both  of  which  participate  in  the 
ozone  destruction  cycle.  Thus,  the  effect  of  reaction  (2.2-28) 
is  to  reduce  both  the  catalytic  role  of  C1X  and  N0x  in  the 
ozone  destruction  processes. 

In  summary,  reactions  (2.2-22),  (2.2-23),  (2.2-24), 

(2.2-25),  (2.2-27)  and  (2.2-28)  provide  the  principal  coupl¬ 

ing  between  the  HOx,  N0x  and  C1X  cycles  in  the  stratosphere. 
Because  of  the  strong  coupling  between  various  chemical 
cycles,  it  is  difficult  to  study  each  cycle  in  isolation. 

It  is  clear,  however,  that  a  thorough  understanding  of  the 
coupling  mechanisms  is  needed  for  the  interpretation  of  field 
data  and  for  model  prediction  of  the  effect  of  the  potential 
perturbations . 

Parallel  to  the  development  of  the  chemical  interaction, 
the  importance  of  dynamic  transport  also  becomes  more  evident. 
As  observation  data  accumulates,  it  becomes  clear  that  the 
seasonal  and  latitudinal  behavior  of  many  gases  cannot  be 
explained  by  photochemistry  alone.  In  the  next  section  we 
will  review  the  treatment  of  dynamic  transport  in  the  dif¬ 
fusion  equation. 


2.3  Treatment  of  Dynamics 


Equation  (2.1-1)  shows  that  the  divergence  of  the  flux 
also  contributes  to  the  time  rate  of  change  of  local  concen¬ 
tration.  In  the  full  3-dimensional  treatment,  <K  is  given 
by  the  standard  form  of  the  Euler ian  fluid  description. 


4> .  =  n .  v 


(2.3-1) 


where  v  is  the  velocity  wind  field  describing  the  general 
circulation.  It  is  implicit  in  (2.3-1)  that  all  trace  species 
follow  the  same  bulk  air  motion. 

Because  of  its  complexity,  a  full  3-D  model  calculation 
places  enormous  demand  on  computer  core  memory  and  computa¬ 
tion  time  as  the  velocity  wind  field  is  a  function  of  space 
and  time.  Somehow,  there  has  to  be  a  trade  off  between 
treatment  of  dynamics  and  chemistry.  In  the  current  3-dimen¬ 
sional  general  circulation  model,  the  chemical  scheme  is  not 
adequate  for  study  of  trace  gas  concentration.  An  alternate 
approach  is  needed  if  one  desires  to  include  complete  chem¬ 
istry. 

In  a  1-D  model,  one  opts  for  a  more  complete  treatment 
of  chemistry  at  the  expense  of  dynamic  description.  It  is 
difficult  to  formulate  the  1-D  model  in  terms  of  global  aver¬ 
aging  from  the  3-D  equation.  Often,  it  is  easier  to  visualize 
the  1-D  model  as  describing  a  one-dimensional  atmosphere.  The 
continuity  equation  for  a  trace  gas  takes  the  form 
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(2.3-1) 


t- 


where  z  is  the  geometrical  altitude,  <t>  is  formally  Wn^  , 
dz 

with  W  =  the  vertical  velocity.  The  flex  can  be  I  (u>u*|IU 
of  as  being  comprised  of  two  parts,  one  due  to  large  scale 
motion  and  the  other  due  to  turbulent  diffusive  motion. 

The  one-dimensional  continuity  equation  implies  that  the 
large  scale  motion  does  not  contribute  to  the  flux.  The 
flux  due  to  turbulent  diffusive  motion  can  be  parameterized 
by  the  eddy  diffusion  coefficients  and  the  gradient  of  the 
mixing  ratio  of  the  trace  gas,  i.e., 

3f  ■ 

=  -K(z)N(z)g—  (2.3-2) 


where  K(z)  is  the  positive  definite  eddy  diffusion  cooffi- 

n. 

cient,  N  (z)  is  air  number  density,  f.  -  -1,  is  the  volume 

1  N 

mixing  ratio.  Note  that  the  flux  <fc  flows  in  opposite  direc¬ 
tion  to  the  mixing  ratio  gradient.  In  the  absence  of  photo¬ 
chemical  interactions,  the  steady  state  solution  is  yiven  by 
f i  us  a  constant  in  altitude.  Equation  (2.3-1)  can  be  re¬ 
written  in  the  form 


i 


at 


3f . 

<-KN— L) 


(2.3-3) 


In  order  to  solve  (2.3-3),  one  needs  a  chemical  scheme 
to  set  up  Qif  and  prescribed  values  for  N(z)  and  K(z).  N(z) 
is  usually  obtained  from  tabulated  values  compiled  from  obser- 
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vations  (U.S.  Standard  Atmosphere1  Supplement  ,  1966).  K(z) 

is  parameterized  by  a  semi-empirical  approach.  Equation 
(2.3-3)  is  solved  wilh  parameterized  values  of  K(z)  for 
species  for  which  the  chemistry  is  well  understood.  K(z)  is 
then  adjusted  so  that  the  computed  altitude  profiles  compare 
favorably  with  observations.  Initially,  methane  has  been 
used  in  this  kind  of  calibration  study  (Wofsy  and  McElroy, 

1973;  Hunten,  1975)  because  of  its  simple  chemistry.  The 
resulting  K(z)  is  then  used  to  solve  for  the  altitude  pro¬ 
files  of  other  trace  gases.  The  computed  results  compare 
favorably  with  observation,  thus  giving  support  to  the  idea 
that  the  same  eddy  diffusion  coefficient  can  be  used  for  all 
trace  gases.  The  eddy  diffusion  coefficient  used  in  the 
present  model  is  taken  from  Wofsy  (1976). 

The  simulation  of  transport  by  vertical  mixing  in  1-D 
models  has  been  applied  with  success  to  gases  like  CH^,  ^O, 

F- 11  and  F-12  for  which  the  stratospheric  concentrations  are 
maintained  by  upward  diffuson  of  ground  level  emissions  to 
balance  the  stratospheric  chemical  removals.  However,  it  is 
found  that  the  1-D  model  has  trouble  accounting  for  all  the 
observed  latitudinal  and  seasonal  behavior  of  the  trace  gas. 
This  is  particularly  true  in  the  case  of  ozone  where  the  oc¬ 
currence  of  a  spring  time  maximum  in  ozone  column  (cf.  Diitsch, 
1971)  is  completely  contrary  to  one’s  expectation  based  on  a 
photochemical  argument. 

Because  of  its  one-dimensional  nature,  the  1-D  model  can¬ 
not  incorporate  horizontal  transports.  However,  the  presence 
of  meridional  transport  is  evident  even  in  the  study  of  the 
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classical  structure  of  the  Hadley  Coll  circulation.  Recent 
measurements  of  CH4  and  N20  reported  by  Ehhalt  (1980)  indicate 
that  1-D  calculations  typically  overpredict  the  stratospheric- 
concentration  of  the  gases  around  the  equator  while  under¬ 
predicting  it  at  high  latitude.  This  is  consistent  with 
the  idea  of  the  upwelling  motion  in  the  equatorial  region 
followed  by  poleward  and  downward  motion  at  higher  latitudes. 

Another  aspect  of  horizontal  transport  comes  from  the 
zonal  winds.  However,  they  are  of  more  interest  in  meteor¬ 
ology  rather  than  photochemistry  since  the  longitudinal 
behavior  of  most  species  are  more  likely  to  be  dominated  by 
the  day-night  variations.  Though  some  questions  remains  as 
to  whether  planetary  waves  could  give  rise  to  net  meridonal 
transport  (cf .  Hartman  and  Garcia,  1979) ,  at  this  stage  of 
the  development,  it  seems  that  the  2-D  model  is  a  good  com¬ 
promise  for  treatment  of  dynamics  and  chemistry. 


2 . 4  Concluding  Remarks 

Within  the  past  five  years,  because  of  the  ever-increas¬ 
ing  computational  power  of  high  speed  computers,  it  has  be¬ 
come  feasible  to  incorporate  a  complete  chemical  scheme  into 
a  2-D  model.  However,  because  of  the  intrinsic  time-depen¬ 
dence  of  the  circulation,  it  is  necessary  to  solve  the  dif¬ 
fusion  equations  via  a  direct  time  integration.  in  contrast, 
the  1-D  formulation  can  be  modified  to  solve  for  the  steady 
state  solution.  Because  of  its  flexibility,  the  1-D  model 
will  remain  useful  as  a  tool  for  identifying  problems  and 
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for  obtaining  first  order  solutions  in  sensitivity  analysis 
In  particular,  the  1-D  model  provides  the  appropriate  test¬ 
ing  ground  for  analyzing  the  ever  more  complex  coupling 
among  the  chemical  species.  On  the  other  hand,  as  more 
global  observational  data  accumulates,  particularly  those 
derived  using  satellite  observations,  it  would  be  increas- 
ingly  important  to  use  a  2-D  model  for  comparisons.  Thus, 
while  it  is  imperative  to  develop  a  more  efficient  2-D 
model,  the  1-D  model  remains  indispensible  as  a  diagnostic 
tool  to  complement  the  2-D  model. 
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3.  the  aer  one-dimensional  model 


3 . 1  Derivation  of  Basic  Equation 

In  this  section,  we  will  outline  the  basic  assumptions 
leading  to  the  derivation  of  the  equations  used  in  the  AER 
1-D  photochemical  diffusive  model  models.  The  structure  of 
the  models,  with  the  diurnal  option  developed  under  the 
present  contract,  will  be  discussed  in  Secton  3.2.  In  addi¬ 
tion,  the  effect  of  diurnal  averaging  on  the  diffusive  species 
calculations  will  be  examined  in  Section  3.3. 

The  basic  equations  in  the  1-D  model  are  the  system  of 
continuity  equations 


(P.  -  L.n. ) } 

l  li 


(3.1-1) 


where  i  is  the  species  index,  and  L^n^  are  the  production 
and  loss  rates  respectively.  The  systems  of  equations  are 
coupled  to  each  other  through  the  production  and  loss  rates. 
Thus,  formally,  the  systems  of  equations  have  to  be  solved 
simultaneously . 

Using  scale  analysis,  one  can  define  the  diffusive  life¬ 
time  tq  and  chemical  lifetime  Tq  with 


(3.1-2) 


Note  that  implicit  in  (3.1-1)  is  the  assumption  that  one 
set  of  K (z)  is  applicable  to  all  species.  Thus,  all  species 
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,  . 


depends  on  the  chemical  reactivity 


have  the  same  i  while 

yi 

of  the  individual  species. 

Using  the  atmospheric  scale  height 

H  =  —  %  6  km 

mg 

as  the  typical  length  scale,  one  estimates  that 

t_  ^  (^-)  1  'v  1  month. 

U  H 

For  atmospheric  species,  t  ranges  from  seconds  to  hundreds 

yi 

of  years.  Thus,  it  is  difficult  to  find  a  single  numerical 
method  that  can  handle  such  a  large  range  of  timescales. 

To  solve  this  problem,  one  makes  use  of  the  chemical 
properties  of  the  atmospheric  species.  The  diffusive  time 
scale  t  provides  a  convenient  standard  to  separate  the  atmo¬ 
spheric  species  into  two  groups.  The  first  group  will  be 
referred  to  as  the  photochemical  species.  This  corresponds 
to  species  having  chemical  lifetimes,  Tq  <<  ;•  (about  1  day 
or  les^  so  that  they  are  always  in  photochemical  equilibrium. 
As  a  result,  the  diffusive  term  can  be  ignored  for  these 
species.  The  rest  of  the  species  will  be  referred  to  as  dif¬ 
fusive  species.  Next  we  examine  the  coupling  among  the  vari¬ 
ous  species.  Because  of  their  short  chemical  lifetimes,  the 
photochemical  species  are  strongly  coupled  to  each  other 
while  only  weakly  coupled  to  the  diffusive  species.  In  con¬ 
trast,  the  diffusive  species  are  only  weakly  coupled  to  each 
other  and  the  photochemical  species.  The  only  exception  may 
be  the  case  of  whose  distribution  directly  affects  trans¬ 
mitted  solar  fluxes  and  modulates  the  local  photolysis  rates. 
Thus,  the  chemical  properties  of  the  atmospheric  species 
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suggest  that  the  systems  of  equations  should  be  treated  as 
two  separate  groups: 


for  photochemical  species 


3n. 

_ i 

at 


[P.  -  L.n. ] 
i  11 


(3.1-3) 


and  for  diffusive  species 


3f.  ,  a  3  f . 

1  =  k{Jt  <KN3ti)  +  (Pi 


at 


at 


Lini)} 


(3.1-4) 


The  partitions  of  the  atmospheric  species  into  the  photochem¬ 
ical  and  diffusive  species  in  the  present  model  are  given  in 
Appendix  A.  The  systems  of  equations  are  to  be  solved  by  an 
iterative  scheme.  Within  each  iteration,  equations  (3.1-3) 
are  to  be  solved  simultaneously  while  (3.1-4)  are  to  be  solved 
one  at  a  time.  This  approach  is  justified  by  the  character¬ 
istic  of  the  chemical  couplings  discussed  before.  The  feed¬ 
back  between  iterations  is  achieved  via  the  revised  produc¬ 
tion  and  loss  rates  arising  from  newly  calculated  species 
densities.  The  algorithm  is  schematically  outlined  in  Fig.  3-1. 

One  major  simplification  arises  in  that  equations  (3.1-3) 
are  uncoupled  in  their  altitude  dependence.  Thus,  they  can 
be  solved  one  altitude  at  a  time.  Equations  (3.1-4)  are 
second  order  in  the  spatial  variables  and  will  require  boun¬ 
dary  values  for  unique  solutions.  In  modeling,  this  usually 
is  done  by  giving  the  flux  at  the  top  of  the  atmosphere  and 
the  value  f ^  or  the  flux  at  the  ground  level.  The  exact 
methods  used  for  solving  (3.1-3)  and  (3.1-4)  will  be  presented 
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in  Appendix  A.  In  the  next  section,  the  time  dependence  oi 
the  equations  will  be  analyzed. 

3. 2  The  Time  Dependence 

Solar  radiation  is  the  ultimate  driving  force  for  atmo¬ 
spheric  chemistry.  As  a  result,  all  the  photochemical  species 
inherit  the  24  hour  periodicity.  For  the  diffusive  species, 
their  lifetimes  are  sufficiently  long  that  they  do  not  re¬ 
spond  to  the  diurnal  variation.  However,  additional  time 
dependence  may  arise  for  boundary  conditions,  e.g.,  the  time- 
dependent  atmospheric  injection  of  F-ll  and  F-12.  In  this 
section,  we  will  concentrate  on  the  treatment  of  the  diurnal 
time  dependence  only. 

If  one  is  not  interested  in  the  diurnal  behavior  of  the 
species,  one  can  obtain  the  diurnal- average  version  of  equa¬ 
tions  (3.1-3)  and  (3.1-4)  by  averaging  them  over  a  period  of 
one  day.  One  then  has 


0 


P . 
x 


L .  n  ■ 


x  x 


(3.2-1) 
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(3.2-2) 


In  the  above  equations  denotes  diurnal  average,  i.e., 

given  h(t,z). 


h(z)  =  jTh?  /  h(t'2,dt  . 
over 
24 

hours 
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In  deriving  (3.2-1)  and  (3.2-2),  it  is  assumed  that  K,  N  are 
not  functions  of  time  and  we  have  ignored  other  time  depen¬ 
dence  for  the  purpose  of  this  discussion.  Next,  we  would 

like  to  point  out  that  the  terms  P.  and  L.n.  consist  of  sums 

111 

of  terms  each  of  which  have  either  one  of  the  following  three 
forms 


J.n.,  k..  n.n.  ,  k..,n.n.n. 

11  13  13  13k  13k 

corresponding  to  photolysis  reaction. 


two- body  and  three- 


body  reaction.  Thus,  one  must  find  seme  way  of  evaluating 


J.n.'  k..n.n.,  and  k..,n.n.n, 
11  13  1  3'  13k  13k 

terns  of  equations. 

Given  g(t,z)  and  h(t,z), 
g*  (t,z)  =  g (t , z) 

h*  (t , z)  =  h(t, z) 


before  one  can  solve  the  sys- 

one  can  define 
g(z) 
h  (z)  . 


Then,  by  definition 

g’  (z)  =  h'  (z)  =  0 

and  gh  =  gh  +  g' h' 


(3.2-3) 


Thus,  the  mean  of  a  product  is  not  equal  to  the  product  of 
the  mean.  The  correction  term  g' h'  cannot  be  estimated  ex¬ 
cept  by  performing  the  full  diurnal  calculation  in  obtaining 
g(t,z)  and  h(t,z).  in  the  diurnal-average  version  of  the 
model,  we  simply  assume  that  the  correction  terms  are  small 
compared  to  the  product.  The  validity  of  this  assumption 
has  been  questioned  (NASA,  1977) .  We  will  discuss  how  good 
this  assumption  is  in  the  next  section. 

In  the  diurnal  version  of  the  model,  equations  (3.1-1) 
are  retained.  The  systems  of  equations 
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(3.2-4) 


9n. (t) 

— at  '  =  pi<t)  -  L.(t)ni(t) 

is  a  non-linear  system  of  differential  equations  in  time. 
Initial  conditions  must  be  specified  before  a  unique  solution 
can  be  obtained.  However,  since  the  initial  conditions  are 
not  exactly  known  for  most  calculations,  we  will  instead 
impose  a  periodic  boundary  condition  on  the  system  by  re¬ 
quiring 


n^(t  +  24  hour)  =  n^(t) 


(3.2-5) 


for  all  species  and  at  all  altitudes.  (Note  that  the  z  de¬ 
pendence  has  not  been  explicitly  included  in  (3.2-4)  since 
the  equations  are  uncoupled  in  altitude.)  The  systems  of 
equations  are  solved  by  first  deriving  a  finite  difference 
scheme  to  propagate  the  solution  ni  (t)  in  time.  The  equa¬ 
tions  are  then  propagated  until  all  species  satisfy  the 
periodic  boundary  conditions.  Further  details  are  presented 
in  Appendix  A. 

The  diurnal  average  version  of  (3.1-2)  is  used  for  the 

diffusive  species  even  in  the  diurnal  option.  The  rationale 

here  being  that  the  concentrations  of  the  diffusive  species 

are  not  going  to  change  appreciably  in  a  day’s  time.  However, 

in  contrast  to  the  diurnal-average  option,  the  term  P  -  L  n 

i  i 

is  computed  exactly  using  the  time-dependent  photochemical 
species  concentrations. 
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3.3  The  Diurnal  Model 


In  this  section,  we  will  discuss  some  of  the  results 
from  the  diurnal  model.  A  detailed  description  oi  the  model 
is  given  in  Appendix  A. 

Clearly,  the  diurnal  model  is  useful  in  interpreting 
the  available  radical  observations.  These  measurements  are 
made  at  definite  local  times  with  definite  sun  angles.  Since 
concentrations  of  radical  species  can  change  by  orders  of 
magnitude  as  a  function  of  the  time,  a  diurnal  calculation 
is  necessary  for  a  meaningful  comparison. 

In  the  present  model,  a  time  step  of  1  hour  is  employed 
during  the  daylight  hours,  2  hours  for  night  time  and  S  hour 
intervals  at  dawn  and  dusk.  Since  an  implicit  scheme  is 
used,  the  choice  of  the  time  step  can  be  made  without  worry¬ 
ing  about  problems  of  numerical  stability.  The  present 
choice  of  time  step  is  adequate  for  the  purpose  of  obtaining 
the  diurnal  behavior  of  the  photochemical  species.  However, 
if  one  is  interested  in  behavior  of  the  species  during  dusk 
and  dawn  (e.g.,  in  case  of  interpreting  limb-scan  data),  the 
time  resolution  during  dawn  and  dusk  should  be  refined. 

Our  calculation  shows  that  many  of  the  radical  species 
have  rather  strong  diurnal  variations.  Fig.  3-2  through 
Fig.  3-4  show  the  diurnal  behavior  of  0,  CIO  and  N02  for 
selected  altitudes.  The  question  arises  as  to  how  the  diurnal 
average  version  of  the  model  approximates  the  actual  behavior 
of  the  radical  concentrations. 

If  ns(z)  is  the  solution  to  the  diurnal  average  equation 
(3.2-1),  one  can  compare  that  to  n (t, z)  where 
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Log  (number  density,  cm'3) 


Figure  3-3 

nC10^t,z*  for  SeJ-ected  Altitudes 
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Figure  3-4 


nN02 (t' z) 


for  Selected  Altitudes 
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n  (  lVzT  =  2  4  ‘iour  /n(L'z)dt 

over 

24 

hours 

with  n(t,z)  being  the  solution  to  (3.2-3).  Fig.  3-5  -  Fig. 

3-7  show  the  comparison  for  0,  CIO  and  NO2  respectively. 

We  note  that  n  (z)  agrees  better  with  n(t,z)  for  the  day 
s 

time  species  than  the  night  time  species.  In  particular, 
ns(z)  grossly  underestimates  n(t,z)  for  night  time  species. 
It  is  clear  why  this  should  be  the  case  as  the  steady  state 
model  represents  the  case  with  continuous  sunlight  and  thus 
would  produce  results  that  approxiamte  the  day  time  number 
densities  for  all  species. 

In  section  3.2,  we  discussed  the  question  concerning 

the  accuracy  of  the  approximation, 

n .  n  .  i  n  .  n  . 
il  il 

In  the  case  of  O^,  one  of  the  terms  in  the  loss  rate  is  pro¬ 
portional  to  the  production  of  nQ(t,z)  and  ncl0(t,z).  With 
proper  diurnal  averaging,  this  term  should  be  given  by 
A(z)  =  nQ(t,z)  nclQ  <t , z)  . 

If  one  ignores  the  correction  term  and  approximates  the  mean 
of  the  product  by  the  product  of  the  mean,  one  gets 
B ( z )  =  nQ (t , z)  ncl0 (t,  zT 

Finally,  in  the  diurnal-average  version,  one  simply  has 

C(z)  =  n  (z)  n  (z) 

0  SC10 

fig.  3-8  shows  a  plot  of  the  three  quantities  as  func¬ 
tions  of  altitudes.  It  is  evident  from  the  curves  that  the 
discrepancy  among  the  various  averaging  methods  could  be  as 
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Figure  3-5 

Plot  of  n  (z)  vs  n(t,z)  for  0 
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Figure  3-6 

Plot  of  n  (z)  and  n(t,z)  for  CIO 
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high  as  factor  of  3.  However,  it  turns  out  that  the  majority 
of  the  terms  in  the  production  and  loss  rates  for  the  diffu¬ 
sive  species  are  not  comprised  of  products  of  photochemical 
species.  Thus,  in  spite  of  the  large  discrepancy  in  the  in¬ 
dividual  terms,  it  only  contributes  to  about  10-20%  of  the 
net  production  and  loss  rates.  As  a  result,  for  the  case  of 
,  the  differences  between  the  concentrations  calculated 
using  the  diurnal-average  version  and  the  diurnal  version 
are  typically  about  10-20%. 
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4. 


COMPARISON  BETWEEN  THEORY  AND  OBSERVATION 


The  comparison  between  model  calculations  with  obser¬ 
vational  data  permits  an  important  test  for  the  validity  of 
the  theoretical  models  used  in  the  assessment  of  any  environ 
mental  problem  including  the  impact  of  fluorocarbons  on 
atmospheric  ozone.  It  may  also  serve  as  a  basis  for  iden¬ 
tifying  any  important  missing  elements  in  our  current  under¬ 
standing  of  the  atmospheric  chemical  and  dynamic  processes, 
should  there  be  serious  discrepancy  between  theory  and  obser 
vations . 

One  of  the  major  difficulties  in  drawing  conclusions 
from  comparing  theory  with  observations  lies  in  the  fact 
that  most  trace  gases  exhibit  large  local  variabilities 
much  of  which  are  not  well  understood.  Although  current 
models  can  simulate  diurnal  variations  of  many  trace  gases, 
they  cannot  account  for  the  local  variability  attributed  to 
dynamics  or  other  atmospheric  fluctuations.  Nevertheless, 
within  the  framework  of  the  1-D  model,  a  number  of  observa¬ 
tions,  notably  the  C10/03,  HN03/N02  and  HF/HC1  data,  cannot 
be  made  entirely  consistent  with  current  theory. 

The  inconsistency  could  simply  be  due  to  the  uncertain¬ 
ties  in  the  currently  adopted  rate  constants  or  the  eddy 
diffusion  profile.  It  may  also  suggest  serious  omission  of 
chemistry  or  dynamics  in  the  models.  In  either  case,  this 
could  have  important  implications  for  the  assessment  of  the 
stratospheric  pollution  problems. 


3 

4.1  Odd  Oxygen  (O (  P)  ,  O3) 

Figure  (4.1-1)  shows  the  calculated  0(3P)  profiles  along 

3 

with  Anderson's  data.  Although  the  theoretical  0(  P)  profile 
lies  within  the  envelope  defined  by  the  data,  it  should  be 

3 

noted  that  the  spread  in  0 (  P)  data  is  relatively  wide 
(-factor  ±2)  in  the  region  below  ~33  km.  Jn  the  context  of 
present  photochemical  theory,  0(  P)  should  be  in  photochemical 
equilibrium  with  O^.  Thus,  for  a  given  solar  zenith  angle, 
the  variations  in  0(3P)  mainly  reflect  the  variations  in 
density.  Quantification  of  variations  are  therefore  impor¬ 
tant  for  ascertaining  the  significance  in  the  observed  varia¬ 
bility  of  O ( 3P) . 

3 

Table  (4.  bl)  shows  the  calculated  0(  P)/0  }  ratio  along 
with  Anderson's  observations  taken  on  12/2/77.  The  calcu¬ 
lated  ratio  is  relatively  insensitive  to  the  diffusion  co¬ 
efficients  and  to  the  photochemistry  of  NO  ,,  HO  and  C1X. 

To  an  excellent  approximation,  the  0.0 ^  ratio  is  given  by 


[0(3P)  ]  =  _ J_3_ _ 

103]  k1L02)(M] 


(4.1-1) 


where  denotes  the  photolysis  rate  of  0^  and  k^  the  rate 
constant  for  the  reaction  0  +  02  +  M  0^  +  M.  Equation  (4.1-1) 
defines  one  of  the  most  basic  relations  in  stratospheric 
chemistry.  While  there  seems  to  be  good  agreement  between 
the  calculated  and  observed  0/0^  ratios,  a  definitive  conclu¬ 
sion  cannot  yet  be  drawn.  It  should  perhaps  be  pointed  out 
that  the  calculated  assumes  some  kind  of  planetary  albedo, 
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TABLE  4.1-1 


OBSERVED  VS.  CALCULATED  10  (  i>  >  J  /  L  O  -  j 

12/2/77  x  =  50° 


ALTITUDE 

OBSERVED 

CALCULATED 

CALCULATED 

(km) 

[0(3P)  ]/  (o3) 

Logan  et  al. 
(1978) 

Sze  et  al. 
(1978) 

4.8 

(-3) 

3.8 

(-3) 

4.5 

(-3) 

3.0 

(-3) 

2.4 

(-3) 

2 . 85 

(-3) 

1.6 

(-3) 

1.6 

(-3) 

1. 8 

(-3) 

1.0 

(-3) 

1.0 

(-3) 

1.1 

(-3) 

6.2 

(-4) 

6.2 

(-4) 

7.17 

(-4) 

3.7 

(-4) 

4.1 

(-4) 

4.5 

(-4) 

2 . 2 

(-4) 

2.6 

(~4) 

2. 9 

(-4) 

1.4 

(-4) 

1.7 

(-4) 

1.89 

(-4) 

9.6 

(-5) 

1.1 

(-4) 

1.05 

(-4) 

7 .7 

(-3) 

7 . 6 

(~  3 ) 

8.04 

(-5) 

3 . 4 

(-3) 

3  .  I 

(-3) 

3 . 3 

(-3) 

4.  3 

(-5) 

3.  4 

(-5) 

3.5 

(-5) 

3.1 

(-5) 

2 . 2 

(-5) 

2.4 

(-5) 

4.8  (-3)  denotes  4.8  x  10 


while  Anderson's  experiment  drd  nut  include  a  measurement 
of  albedo.  Thus  the  agreement  could  be  fortuitous. 

Figure  (4.1-2)  shows  the  calculated  O^  protiles  along 
with  observations.  The  behavior  below  26  km  is  very  sen¬ 
sitive  to  HO^  chemistry.  In  our  earlier  calculation,  the 
calculated  concentrations  below  26  km  as  well  as  the  total 
column  density  are  20%  higher  than  observation  as  a  result  of 
the  revised  rate  constant  (Howard  and  Evenson,  1977)  for  the 
reaction. 


H02  +  NO  -*■  OH  +  N02  .  (4.1-2) 

Recent  revision  in  another  rate  constant  for  the  important 
reaction  of  H02  with  0^  removes  most  of  the  discrepancy 
between  the  calculated  and  the  observed  ozone  abundances 
below  26  km. 

The  calculated  O^  concentration  above  40  km  however, 
still  seems  to  lie  somewhat  below  most  observation  data 
particularly  those  of  Riegler  et  al.  (1977).  It  is  unfor¬ 
tunate  that  the  apparent  discrepancy  between  various  sets  of 
0  ^  data  (Kruger  and  Minzner,  1976;  Wantanabe  and  Tohmat.su,  1976; 
Hering  and  Borden,  1967;  Riegler  et  al.,  1977)  are  not  yet 
resolved.  It  would  seem  that  simultaneous  measurements  of 

3 

0  (  P)  ,  and  Oil  would  permit  a  relatively  straight  forward 
test  for  the  odd  oxygen  balance  above  44  km  (Sze,  1979),  a 
region  where  removal  of  odd  oxygen  is  dominated  by  the 
reactions  0  +  OH  •+■  H  +  02<-  0  +  H02  -*•  OH  +  02  and  0  +  03  ■+  20 2« 


.j  <w  • 
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ALTITUDE  (km) 


4.2  Ilydrocjon  Radicals  (Oil,  110^,) 


Figures  (4.2-1)  and  (4.2-2)  show  the  calculated  OH  and  HO^ 
profiles  at  solar  zenith  angle  Y  =  40°  along  with  Anderson's 
data.  Agreement  between  theoretical  and  observed  |noi lies 
seems  to  be  reasonably  good.  Data  for  Oil  and  llo^  below  30 
km,  however,  are  still  lacking.  The  CI1  concentration  in 
this  region  (z  <  30  km)  may  only  be  inferred  indirectly  from 
simultaneous  measurements  of  certain  ratios  (e.g.  HF:HC1; 
hNO ^ :  NO2 ;  C.\0:HC1)  or  from  certain  atmospheric  species 
(e.g.  C2^^)  whose  altitude  distributions  are  extremely 
sensitive  to  OH  abundance. 

4. 3  Free  Chlorine  Species 

Figure  (4>l)shows  the  calculated  CIO  profiles  along  with 
observations.  There  are  six  CIO  profiles  reported  by  Anderson 
et  al.  (1979)  and  one  by  Menzies  (1979).  An  excellent  discus¬ 
sion  on  the  CIO  data  and  the  implications  for  stratospheric 
chemistry  is  given  by  Anderson  et  al.  (1980) .  We  shall  sum¬ 
marize  here  some  of  the  interesting  features  regarding  the 
CIO  distribution. 

Four  of  Anderson's  CIO  profiles  measured  in  Call  and 
winter  lie  below  1  ppbv,  whereas  two  of  the  summer  profiles 
exceed  1  ppbv  above  30  km.  During  the  July  14,  1977  experi¬ 
ment,  simultaneous  OH  and  0^  concentrations  were  measured 
along  with  CIO.  The  observed  high  CIO  data  (with  peak  mixing 
ratio  '8  ppbv)  and  normal  0^  are  not  explicable  in  the  present 
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OH  MIXING  RATIO 


Figure  4.2-1 

Calculated  altitude  profiles  for  OH  mixing  ratio  corresponding 
to  41°  zenith  angle.  The  data  from  Anderson  (see  NASA,  1979) 
are  included  for  comparison.  The  observations  were  made  at 
32  °N.  The  zenith  angles  during  observations  are  given  in  the 
figure.  Note  that  the  observations  with  zenith  angle  of  80° 
have  been  adjusted  by  factor  of  two  for  comparison  purposes. 
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ALTITUDE  (km) 


ho2  mixing  ratio 


Figure  4.1-2 

Calculated  altitude  profiles  lor  mixiny  ratio  corresponding 

to  41°  zenith  angle.  The  data  from  Anderson  (see  NASA,  1979) 
and  Mihelcic  et  al.  (1978)  are  included  for  comparison.  The 
observations  were  made  at  30°N  and  53°N  respectively.  The 
zenith  angles  during  observations  are  given  in  the  figure. 


4-9 


CK)  VOLUME  MIXING  RATIO 


Figure  4.3-1 


Calculated  altitude  profiles  of  CIO  mixing  ratio  corresponding 
41°  zenith  angle.  The  observations  from  Anderson  et  al.  (1980) 
and  Menzics  (1979)  are  included  for  comparison.  Zenith  angles 

for  the  observations  are 
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theoretical  framework.  Anderson  et  al.  (1980)  argued  that 
the  large  values  of  CIO  could  not  be  attributed  to  experi¬ 
mental  error.  They  concluded  that  the  observations  implied 
either  (a)  that  significant  quantities  of  CIO  were  ini'-'  t.-u 
into  the  study  region  prior  to  their  measurement  ,  on  a  time- 
short  compared  with  the  chemical  response  t  nue  for  odd 
oxygen  (about  a  day  at  40  km),  or  (b)  that  the  reaction  ot 
CIO  with  0  was  not  rate  limiting  in  Liu.'  closure  of  tin 
chlorine  catalytic  cycle,  or  (c)  that  there  exist  sources 
for  odd  oxygen  proportional  to  CIO^  which  can  compete 
directly  with  photolysis  of  in  the  middle  and  upper 
stratosphere . 

Another  serious  discrepancy  revealed  by  Anderson  et 
al.'s  data  concerns  the  rapid  fall  off  of  CIO  mixing  ratio 
below  30  km.  One  straight  forward  way  to  reconcile  the 
lower  stratospheric  CIO  data  is  to  suppress  the  Oil  concen¬ 
tration  below  30  km. 

Figure  (2.6)  shows  the  calculated  C10N02  profile  along 
with  recent  data  by  Murcray  et  al.  (1979).  The  apparently 
good  agreement  between  calculated  and  observed  CluNO^,  concen¬ 
trations  may  be  fortuitous  in  view  of  the  paucity  of  C  IONo^, 
data  and  uncertainties  in  CIONO^  chemistry.  For  instance, 
there  is  no  direct  evidence  that  the  reaction  CIO  with  No^ 
yields  CIONO2  or  CINO^  isomer  as  the  primary  product. 
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4.4  HF : HC1  and  11N03:NC>2  Ratios  as  Indicators  for 
Stratospheric  OH 

Fiyures  (4. 4-  1)  and  (4  4-2)  show  the  calculated  lU':llCi  and 
HN03:N02  ratios.  These  two  ratios  are  quite  sensitive  to 
stratospheric  OH  density  (NASA,  1977;  McConnel  and  tvuns, 
1978;  Sze,  1978,  1979).  Comparison  between  calculated 
ratios  and  observations  permits  an  indirect  check  on  the 
predicted  stratospheric  OH  abundance  below  30  km,  where  OH 
data  are  still  lacking. 

The  calculated  HNO,:NO-,  ratio  appears  to  lie  above  the 

Ihno3) 

observed  ratio  by  about  a  factor  2-4.  The  ratio  — — —  in 


[no2] 


the  model  is  set  mainly  by  a  balance  of 


N02  +  OH  +  M  ->■  HN03  +  M 


hv  +  HNO . 


k2 


OH  +  NO. 


and 


OH  +  HN03  H20  +  N03 


(4.4-1) 
(4.4-2) 
(4  .  4-3) 


with  (44-2)  more  important  than  (44-3)  It  should  be  defined 
to  adequate  precision  in  the  lower  stratosphere  by 


[HN03]  k1(OH][M3  kjlomiM] 

[N02]  =  J3  +  k2[OH]  " 


The  rate  constants  k^  and  J3  are  relatively  well  known. 
It  is  difficult  therefore  to  conceive  of  errors  in  compu¬ 
tation  which  might  cause  ratios  obtained  from  (4  4-4)  to  be 
high  by  a  factor  as  large  as  4.  The  discrepancy  might  be 
taken  to  indicate  a  major  error  in  computed  values  for  OH, 
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ALTITUDE 


U  ~  O 


ALTITUDE  (km) 


hno3/no2  concentration  ratio 


Figure  4.4-2 

Ratio  of  HNO3/NO2  as  computed  from  models  compared  with  data 
The  computed  values  are  for  the  ratios  at  sunset  for  30°N. 

The  data  points  are 

O  Evans  et  al.  (1976)  sunset  59°N 
0  Harris  et  al.  (1976,  1978)  noon  44°N 
□  Fontanella  et  al.  (1975)  Sunset  45-50°N 
A  Loewenstein  et  al.  (1978)  daytime  20-40°N 
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m  such  a  sense  that  model  values  for  Oh  in  the  lower  strato¬ 


sphere  might  be  too  high.  A  reduction  in  OH  would  be  bene¬ 
ficial  also  for  CIO.  Reduction  in  OH  would  tend  to  diminish 
the  source  of  Clx  from 

OH  +  HC1  --  H20  +  Cl  (4.4-5) 

and  if  the  change  in  OH  occurred  primarily  at  low  altitude, 
one  might  expect  to  find  better  agreement  between  model  and 
observed  gradients  of  CIO  as  well  as  the  ratio  of  HF:HC1. 

On  the  other  hand,  the  model  appears  to  give  satisfactory 
results  for  Oil  and  110^  at  higher  altitudes  as  shown  in 
figures  (43-1)  and  (4.2-2). 

There  are  several  ways  that  we  might  overestimate  the 
OH  concentration  below  30  km.  For  example,  considerable 
uncertainty  still  exists  in  the  rate  constant  for  the  reac¬ 
tion 

OH  +  H02  -  H20  +  02  (4.4-6) 

Hubora Lory  studies  by  various  groups  (liochanadel  et  al., 

1972;  burrows  et  al.  ,  1978;  Chang  and  Kaufman,  1978;  DeMore , 

.197  9)  indicate  that  k^.  may  be  pressure-dependent.  Should  k^ 
be  10  cm'*  s  ^  below  28  km  and  •'2x10  ^  cm^  s  ^  above 
14  km,  much  of  the  inconsistency  between  calculated  and  ob¬ 
served  ratios  of  HF/HC1  and  HN03/N02  may  be  resolved,  without 
contradicting  Anderson's  OH  data  above  30  km.  It  is  clear 
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that  the  faster  rate  constants  tor  the  reactions  of  Oil  with 


HN03  and  H02NC>2, 

k7 

OH  +  HN03  H20  +  N03  (4.4- 

k8 

OH  +  H02N02  ->■  H20  +  N02  +  02  (4.4- 

could  also  reduce  the  OH  in  the  lower  stratosphere.  Recent 
measurement  by  Wine  et  al.  ( 1980)  yives  the  following  expres¬ 
sion  for 

k^  =  1.52x10  ^  exp(^T')  .  (4.4- 

Note  that  reaction  (4.4-7)  has  a  neyat Lve  activation  oncryy 

of  about  1.3  kcal  mol  At  stratospheric  temperatures,  the 

new  rate  constant  k^  (Wine  et  al.,  1980)  is  about  a  factor  of 

-14 

3  faster  than  the  temperature  independent  value  of  8.5x10 
cm^  s  recommended  earlier  by  NASA  (1979).  If  Wine  et  al.  's 
results  are  confirmed  and  H20  is  the  major  product  of  reac¬ 
tion  (4.4-7),  then  much  of  the  discrepancy  between  observed 
and  calculated  ratios  of  11F:HC1  and  iiN03:N02  may  be:  icsolved. 

Reaction  (4.4-8)  may  further  reduce  OH  in  the  low  strato- 

-12  3  -1 

sphere  if  kg  exceeds  1x10  cm  s  .  Obviously,  accurate 
data  for  k, ,  k_  and  kQ  are  needed  for  modeling  the  strato- 
spheric  OH  which  plays  a  central  role  in  the  03  problem. 
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SULFUR  STUDY 


5 . 1  Background 

The  interests  in  atmospheric  sulfur  chemistry  originated 
with  the  observation  of  the  stratospheric  sulfate  (Junge) 
layer  (Junge  et  al.,  1961).  Studies  showed  that  the  sulfate 
layer  contributes  to  the  net  albedo  of  the  atmospheric-earth 
system  through  absorption  and  backscattering  of  solar  radia¬ 
tion  (Mitchell,  1971;  Chylek  and  Coakley,  1974;  Pollack  et 
al.,  1976).  Consequently,  any  significant  perturbation  of 
the  atmospheric  sulfur  budget  may  cause  changes  in  the  sul¬ 
fate  layer  with  possible  global  climatic  impact.  Thus,  it 
is  important  to  understand  the  nature  of  the  sulfur  sources 
that  maintain  the  present  sulfate  layer. 

Junge  suggested  that  S02  is  the  gaseous  precursor  of 
the  sulfate  aerosol  and  that  the  stratospheric  S02  concentra¬ 
tion  is  maintained  by  upward  diffusion  from  the  troposphere 
(1974) .  In  this  way,  the  stratospheric  sulfur  budget  is 
related  to  the  tropospheric  budget.  In  the  meantime,  inter¬ 
est  in  the  tropospheric  budget  is  growing  because  of  the 
emerging  problem  in  visibility  reduction,  acid  rain  and  res¬ 
piratory  effects  associated  with  enhanced  SO^  concentration 
due  to  S02  emission  from  anthropogenic  activities. 

With  the  projected  increase  of  coal  usage  in  the  future, 
it  becomes  increasingly  desirable  to  ascertain  the  impact  of 
S02  emission  on  the  global  budget.  A  natural  first  step  in 
the  study  is  the  establishment  of  the  atmospheric  clean  air 
sulfur  budget. 


5 . 2  Natural  Sulfur  Budget 


During  the  past  contract  year,  we  initiated  a  study  of 
the  clean  air  sulfur  budget.  In  this  section,  some  of  the 
major  results  will  be  outlined.  The  reader  is  referred  to 
the  papers  cited  for  more  detailed  discussion. 

The  chemistry  of  the  atmospheric  sulfur  species  and 
available  observations  are  summarized  in  Sze  and  Ko  (1980b) . 
The  measured  tropospheric  concentration  of  COS  of  .5  ppbv 
makes  COS  the  most  abundant  sulfur  species  in  the  atmosphere. 
One  major  conclusion  of  our  study  is  that  reduced  sulfur  com¬ 
pounds  such  as  l^S,  CS2»  COS,  CH^SCH^  may  play  an  important 
role  as  sources  for  atmospheric  S02  and  S04 . 

As  suggested  by  Crutzen  (1976),  photolysis  of  COS  in 
the  stratosphere  could  provide  a  source  for  stratospheric 
S02 .  Our  studies  showed  that  the  reaction  of  CS2  with  O 
could  also  provide  a  source  for  S02  of  comparable  strength. 
Furthermore,  our  model  calculations  showed  that  the  in  situ 
production  of  SO2  from  oxidation  of  COS  and  CS2  can  explain 
why  S02  is  well  mixed  in  the  lower  stratosphere  (cf.  Sze  and 
Ko,  1979b,  1980a). 

Our  studios  also  indicated  that  CS2  may  be  a  precursor 
for  atmospheric  COS  (Sze  and  Ko,  1979a;  Ko  and  Sze,  1980)  . 
Using  the  chemical  scheme  discussed  in  Sze  and  Ko  (1980b), 
we  argued  that  an  atmospheric  input  of  27  MT(S)/year  in  the 
form  of  reduced  sulfur  compounds  can  account  for  the  observed 
SO2  and  S04  tropospheric  concentrations  at  remote  areas. 

This  suggests  that  the  anthropogenic  input  of  65  MT(S)/year 
could  have  a  significant  impact  on  the  global  sulfur  cycle. 
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6 .  FFFKCT  Of  MAN’S  ACTIVITIES  ON  STRATOSPHERIC  O  ^ 

Several  human  activities  which  could  appreciably  per¬ 
turb  the  stratospheric  ozone  layer  have  been  identified  in 
recent  years.  These  include  the  use  of  high-flying  aircraft, 
halocarbons,  and  nitrogen  fertilizers.  Among  the  various 
possible  threats  to  the  ozone  layer,  perhaps  the  release  of 
fluorocarbons  appears  to  be  the  greatest  (NAS,  1980).  This 
section  will  present  the  most  recent  results  concerning  the 
0^  perturbations  associated  with  the  release  of  F-ll  and  F-12 
(Section  6.1),  CH^Br  (Section  6.2)  and  the  injection  of  NOx 
by  high-flying  aircraft  (Section  6.3). 

6.1  0-j  Depletion  due  to  Fluorocarbons 

As  discussed  in  Section  2,  chlorine  radicals  play  an 
important  role  in  stratospheric  ozone  budget.  One  poten¬ 
tially  large  source  of  stratospheric  chlorine  is  the  photol¬ 
ysis  of  CFC13  (F-ll)  and  CF2C12  (F-12), 

CFC13  +  hv  ■+  CFC12  +  Cl  (0.1-1) 

and 

CF2C12  +  hv*  ►  CF2C1  +  Cl  (6.1-2) 

as  identified  by  Molina  and  Rowland  (1974). 

The  two  compounds  (F-ll  and  F-12)  are  believed  to  be 
entirely  man-made.  The  production  rates  of  F-ll  and  F-12 
jl n  the  period  (1950-1979)  are  summarized  in  Table  0.1-1. 
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A  Tabulation  of  the  Total  Production 
of  F-ll  and  F-12  on  an  Annual  Basis 
for  the  Period  1950-1979  (Source  CMA) 


,  F-ll 

F-i2 

Year 

(10°  pounds) 

(10  pound 

1950 

14.6 

76.2 

1951 

20.0 

7  9.9 

1952 

29.9 

82.1 

1953 

38.  1 

102.5 

1954 

46.1 

108.3 

1955 

57.9 

127.0 

1956 

71.6 

151.4 

1957 

74.8 

163.5 

1958 

65.1 

161.9 

1959 

78.4 

193.1 

1960 

109.6 

219.2 

1961 

133.3 

2  3  9.2 

1962 

172.2 

282.4 

196  3 

205.7 

322.8 

1964 

244.9 

375.0 

1965 

270.8 

419.0 

1966 

310.  9 

476.6 

1967 

352.2 

535.2 

1968 

403.7 

589.7 

1969 

479.0 

655.4 

1970 

525.0 

707.9 

1971 

580.2 

753.0 

1972 

676.5 

837.5 

1973 

769.6 

933.3 

1974 

815.1 

976.2 

1975 

692.4 

839.9 

1976 

749.2 

905.5 

1977 

706.5 

844.0 

1978 

680.9 

82  0.3 

1979 

638.2 

787.4 
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These  data  were  recently  compiled  by  the  Chemical  Manufuc- 
turers  Association. 

No  tropospheric  sinks  for  F-ll  and  F-12  have  yet  been 
identified.  Thus,  the  steady  state  concentrations  of  these 
compounds  may  reach  as  high  as  1-2  ppb  from  their  present 
levels  of  .1-.2  ppb  in  the  next  150  years  if  the  present 
production  rates  prevail. 

Several  groups  have  calculated  the  steady  state  ozone 
depletions  due  to  the  release  of  F-ll  and  F-12.  The  results 
are  summarized  in  Table  6.1-2.  The  calculated  reductions  in 
ozone  range  from  15-18  percent  based  on  the  rate  data  recom¬ 
mended  by  NASA  (1979) .  Note  that  the  current  model  predicted 
ozone  depletion  values  are  about  a  factor  of  two  higher  than 
that  published  in  the  NAS  (1976)  report.  The  revisions  in 
the  calculated  ozone  depletion  is  mainly  caused  by  the  up¬ 
ward  revisions  in  the  rate  constants  for  two  key  reactions, 


H02 

+  NO  -*■ 

OH  + 

no2 

(6  . 

1-3) 

H02 

+  °3  " 

OH  + 

202  . 

(6. 

1-4) 

Both  reactions  (6.1-3)  and  (6.1-4)  enhance  the  OH  con¬ 
centrations  below  32  km  with  a  consequent  increase  in  the 
catalytic  role  of  C1X  as  a  sink  for  O^. 

Figure  6.1-1  shows  the  steady  state  locaj[  ozone  reduc¬ 
tions.  It  should  be  emphasized  that  while  the  relative 
change  in  O^  concentration  peaks  at  about  40  km,  the  bulk 
of  ozone  reduction  occurs  below  12  km,  a  region  of  consider¬ 
able  complexity  both  in  terms  of  dynamics  and  chemistry. 
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Table  6.1-2 


Steady-State  Pertubatrons  for  1975  Emissions  Rates 
of  Fluorocarbons  11  and  12* 

(NASA,  1979) 


I 

i 

I 

I 

i 

I 

! 

i 

! 


Modeling 

Group 

AC1X 

ppbv 

-AO  3  (col) 
U) 

-AO  3  (40  km) 
'(%) 

-AO 3 (20  km) 

(  o  ) 

-AO3  (col 

AC1X 
('i»/  ppbv 

bbb 

5.5 

15 

38 

17 

2.7 

UuPoUt 

6 . 2 

lb.  3 

46.  i 

16.5 

3.0 

LaKC 

6.7 

16.2 

43 

h- 

CO 

2 . 4 

KUA 

5.5 

16 

42 

17 

2. 9 

Harvard 

6  .  b 

lb 

41 

19 

2.7 

Cal  Tech 

6.  y 

lb 

50 

19 

3.  1 

WOAA 

5  .  U 

lb 

3y 

17 

3. 1 

GSFC 

7.0 

18 

39 

18 

2.6 

GSFC-UWCJ 

o.  y 

15 . 4 

53 

9 

2 . 2 

t\l* Li 

5  .  U 

lo  .  / 

3  2 

14 

2 . 9 

* 


Fluorocarbon  11  at  9.4  x  10^  and  Fluorocarbon  12  at 

6  *2  - 1 
L3.0  x  10  molecules  cm  sec 
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ALTITUDE 


AO -i  (  in  10"mo  ecules  cm-3) 


AO3  ( in  percentage ) 


figure  0.1-1 

Calculated  steady  state  local  reduction  of  ozone  assuraincj  constant 
10/ /  release  rates  ol  l'-ll,  K~  IZ .  The  dotted  curve  col  responds  lo 
the  reduction  in  number  density  (upper  scale)  while  the  solid  curve 
gives  the  local  percentage  reduction  (lower  scale). 


W;-  ia. 


Note  that  this  is  also  the  region  where  we  found  significant 
discrepancy  between  models  and  observations  (see  Section  4). 
Thus,  the  calculated  changes  in  ozone,  particularly  those 
calculated  below  32  km,  should  be  viewed  with  great  caution. 


6 . 2  Perturbation  due  to  Bromine  Compounds 

Bromine  can  affect  ozone  through  a  number  of  catalytic 
cycles : 


(6.2-1) 

(6.2-2) 

(6.2-3) 

The  sequence  (6.2-1)  is  equivalent  to 

0  +  03  >  202 

while  sequences  (6.2-2)  and  (6.2-3)  are  equivalent  to 

203  -*  30 2  • 


and 


Br  +  03 

-> 

BrO  +  02  1 

BrO  +  0 

Br  +  02  J 

Br  +  03 

->• 

BrO  +  02  \ 

BrO  +  BrO 

-> 

2Br  +  02  f 

Br  +  03 

-> 

BrO  +  02 

Cl  +  O-, 

-> 

CIO  +  o2 

CIO  +  BrO 

-y 

Cl  +  Br  +  O 

The  troposphere  contains  a  variety  of  brominated  organic 
compounds  including  CH^Br ,  CH2Br2,  CHBr  3 ,  CHClBr2  and  C2H2Br2 


6-t> 


Cuncen- 


(Singh  et  al. ,  1977) ,  with  CH^Br  the  must  abundant, 

trations  of  CH^Br  range  from  1  to  300  pptv,  with  values  in 
clean  continental  air  between  5  and  10  pptv.  Tin-  lifetime 
of  CH-jBr,  as  set  by  the  reaction  with  OH,  is  sufficiently 
long  (of  about  3  years)  to  insure  that  the  gas  should  be 
relatively  well  mixed  in  the  troposphere  and  act  as  a  bro¬ 
mine-carrier  to  the  stratosphere. 

Recently  Yung  et  al.  (1980)  argued  that  sequence  (6.2-3) 
may  be  potentially  more  important  because  of  the  possibility 
of  further  increase  in  CIO  concentrations  in  the  stratosphere. 
They  argued  that,  with  8  ppb  of  ClX  (which  corresponds  to  a 
steady  state  stratospheric  ClX  concentration  should  emission 
of  F-ll  and  F-12  continue  at  present  rates),  reaction  se¬ 
quence  (6.2-3)  could  deplete  an  additional  2%  column  ozone 
if  the  background  level  of  BrX  in  the  stratosphere  is  about 
20  ppt.  The  impact  of  bromine  on  ozone,  however,  would  be 
largest  in  the  lower  stratosphere  and  would  involve  nonlinear 
coupling  between  BrX  and  ClX  cycles.  The  calculated  concen¬ 
trations  of  BrO  may  be  subject  ot  larger  errors  since  the 
photolysis  rate  of  BrO  is  not  known  within  a  lector  of  4. 

The  perturbation  results  concerning  the  brominated  compounds 
should  be  considered  therefore  subject  to  the  usual  caveats 
as  noted  in  Section  6.1. 
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b.J  NO  Injection  Studies 
x  J _ 

Earlier  model  calculations  (Crutzen,  1970;  Johnston, 

1971;  McElroy  et  al. ,  1974;  CIAP,  1974)  indicate  that  a 
potential  fleet  of  operations  of  supersonic  aircraft  as 
contemplated  by  the  United  States  in  1970  (500  aircraft 
flying  approximately  7  hours  per  day  at  17-18  km)  could 
lead  to  major  reduction  in  the  stratospheric  column 
abundance  and  thus  cause  an  increase  in  the  flux  of  ultra¬ 
violet  radiation  reaching  the  earth's  surface.  This  could 
result  in  a  variety  of  environmental  consequences,  including 
a  possible  increase  in  the  incidence  of  skin  cancer  (McDonald, 
1971)  . 

Removal  of  by  aircraft  injectant  NO^  radicals  is 
primarily  due  to  the  pair  of  reactions 

NO  +  03  -+  N02  +  02  (6.3-1) 

followed  by, 

N02  +  O  ■>  NO  +  02  (6.3-2) 


The  natural  source  for  NOx  in  the  stratosphere  is  thought  to 
emanate  from  the  reaction  of  0(1D)  with  N20  (Nicolet  and 
Vergison,  1971;  Crutzen,  1971;  McElroy  and  McConnell,  1971), 

N20  +  0(iD)  -►  2NO  .  (6.3-3) 
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Much  of  the  revision  in  model  predicuons  since  1976 
may  be  ascribed  to  a  change  in  the  rate  constant  for  the 
reaction  of  NO  with  H02, 

NO  +  H02  -*■  N02  +  OH  (6.3-4) 

Recent  measurements  (Howard  and  Evenson,  1977)  indicate  that 
the  rate  constant  for  reaction  (6. 3-4 )  is  faster  than  previously 
thought,  by  about  a  factor  of  30.  An  increase  in  the  rate 
constant  for  reaction  (6. 3-4)  tends  to  shift  the  equilibrium 
in  HO^  from  HO 2  toward  OH.  Below  30  km,  removal  of  ozone  by 
HO^  radicals  proceeds  mainly  by, 

OH  +  03  -  llO 2  +  02  (6. 3- 5a) 

H02  +  03  -*  OH  +  20 2  .  (6. 3- 5b) 

We  may  note  that  reactions  (6.3-5a)  and  (6.3-4)  followed  by 
photolysis  of  N02 , 

n°2  +  hv  ->•  NO  +  O  (6.3-6) 

do  not  affect  odd  oxygen  removal.  Thus  addition  of  NO  would 
reduce  the  catalytic  role  of  HOx  as  described  by  <6.3-5a,b). 
Furthermore,  an  increase  in  OH  causes  an  increase  in  the 
rate  for  the  reaction 

OH  +  N02  +  M  f  HN03  +  M  ,  (6.3-7) 
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with  subsequent  decrease  in  the  concentration  of  the  free 


nitrogen  radicals  (NO  +  NO2)  and  reduction  in  the  efficiency 
of  the  nitrogen  cycle  as  a  sink  lor  odd  oxygen. 

In  this  section,  we  present  the  result  of  a  model  calcu¬ 
lation  simulating  the  effect  of  NO^  injection  in  the  strato¬ 
sphere  due  to  aircraft  operation.  The  calculations  are  done 
using  the  diurnal  average  version  of  the  AER  1-D  photochem¬ 
ical  model  with  the  rate  data  as  taken  from  NASA  (1979) .  As 
discussed  in  Section  3  and  Appendix  A,  the  model  atmosphere 
is  taken  from  U.S.  Standard  Atmosphere  Supplement  (1966)  and 
the  eddy  diffusion  coefficients  from  Wofsy  (1976). 

In  the  study,  we  assume  that  all  N0x  are  injected  in 

the  form  of  NO.  This  is  a  valid  approximation  as  the  CIAP 

studies  (CIAP,  1975)  have  shown  that  90-95%  of  the  NO  is 

x 

emitted  in  the  form  of  NO.  It  is  estimated  that  1,000  air¬ 
crafts  flying  7  hours  per  day  will  deposit  1.7  MT  (NOx)  yr_1 
in  the  stratosphere.  Assuming  the  emitted  NOx  to  be  uniformly 
distributed  over  a  1  km  thick  layer,  we  estimate  the  equiv¬ 
alent  injection  rate  of  NO  to  be  1.4xl03  molecules  (NO)  cm-3 
s  .  Two  separate  studies  were  performed  assuming  injection 
at  15  km  and  20  km  respectively. 

Figures  6.3-1  and  6.3-2  show  the  resulting  changes  in 
O^  end  NOx  profiles  as  functions  of  altitude  for  the  two 
studies.  Our  results  also  show  a  net  increase  in  0^  column 
density  of  1.7%  and  4.0%  for  the  15  km  and  20  km  injections 
respectively . 

The  following  points  have  to  be  taken  into  account  in 
interpreting  the  1-D  model  results.  There  is  little  doubt 
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ALTITUDE(km) 


ANOX  ( 10^  molecules  cm-3) 


AO3  ( 1010  molecules  cm-3) 


Figure  6.3-1 

Calculated  changes  in  0  and  NO  profiles  for  model  A 
due  to  an  Injection  of  1. 4x10"  nolcculea  (NO)  cn~2  g-l  a£  15-J6  (cm. 
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that  OH  plays  a  pivotal  role  in  stratospheric  chemistry  and 
that  exact  responses  of  the  distribution  to  the  N0x  injec¬ 
tion.  While  the  kinetic  data  base  for  H0x  reactions  has 
been  significantly  improved  over  the  past  years,  remaining 
uncertainties  in  H0x  chemistry  still,  perhaps,  represent  the 
largest  source  of  error  for  model  predictions.  In  fact, 
a  comparison  between  model  calculations  and  observations 
reveals  several  significant  discrepancies  which  might  be 
attributed  to  errors  in  the  calculated  OH  concentrations  in 
the  altitude  region  15-35  km.  It  was  argued  elsewhere  (Sze, 
1978,  1979)  that  significantly  lower  stratospheric  OH  con¬ 
centrations  than  those  calculated  by  current  models  are 
needed  to  account  for  the  observed  gradients  of  CIO  (Ander¬ 
son  et  al. ,  1979)  and  for  the  observed  ratios  of  HNO^/NO^ 
(NASA,  1977  ,  1979;  McConnell  and  Evans,  1978)  and  liF/llCl 
(Sze,  1978). 

Another  area  of  major  uncertainty  concerns  atmospheric 
transport  of  trace  species.  Current  one- dimensional  models 
parameterize  vertical  transport  by  the  so-called  eddy  diffu¬ 
sion  coefficients  which  were  mainly  derived  from  observation 
of  N 20  and  CH^.  These  models  therefore  ignore  horizontal 
transport,  while  the  natural  distribution  of  ozone  exhibits 
significant  latitudinal  and  seasonal  variations. 

It  is  probably  desirable  to  perform  similar  studies 
using  a  two-dimensional  model  once  some  of  the  uncertainties 
in  the  HOX  chemistry  can  be  removed. 
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7. 


THE  AER  TWO-DIMENSIONAL  MODEL 


7 . 1  Current  Status 

The  development  of  the  AER  2-D  model  started  in  the 
summer  of  1979.  The  approach  to  dynamic  treatment  follows 
that  of  the  Oxford  Model  (Harwood  and  Pyle,  1975,  1977). 

The  model  is  designed  ultimately  to  calculate  all  dynamic 
variables  together  with  atmospheric  species  concentrations 
incorporating  feedback  mechanisms  between  chemistry,  radia¬ 
tion  and  dynamics . 

In  the  early  stage  of  the  development,  we  have  set  up 
the  program  to  solve  the  zonal-mean  diffusion  equation  with 
pre-calculated  dynamic  variables  that  are  specified  extern¬ 
ally  as  functions  of  time  and  space.  The  values  are  taken 
from  earlier  work  of  Harwood  and  Pyle  (1977) .  This  model 
was  used  to  perform  inert  tracer  studies  to  assure  the  proper 
treatment  of  the  equation. 

Our  next  step  was  to  interface  a  photochemical  scheme 
with  the  model.  We  decided  to  build  up  the  chemistry  in 
stages.  The  first  chemical  scheme  corresponds  to  the  clas¬ 
sical  Chapman  scheme.  We  next  added  hydrogen  chemistry  and 
finally  implemented  the  full  °x-H°x-NOx-Clx  chemistry  to 
simulate  the  present  atmosphere.  It  is  hoped  that  by  exam¬ 
ining  the  behavior  of  in  each  case,  one  can  get  a  better 
idea  of  the  role  of  chemistry  and  dynamics  in  determining 
trace  gas  distribution. 

In  the  next  section,  we  will  discuss  the  derivation  of 
the  basic  equations  used  in  our  zonal  model.  Emphasis  will 


be  pul  on  the  physical  assumptions  that  are  adopted  m  order 
to  obtain  a  set  of  manageable  equations  from  the  full  three- 
dimensional  primitive  equations.  A  brief  description  oi 
the  model  can  be  found  in  Appendix  3.  The  results  of  the 
numerical  calculations  will  be  discussed  in  Section  7.1. 

7 . 2  Basic  Equations  for  Zonal-Mean  Model 

The  basic  equations  used  in  zonal  models  can  be  consid¬ 
ered  as  obtained  from  the  full  three-dimensional  equations 
by  zonal  averaging.  In  what  follows,  we  will  begin  with  a 
discussion  of  the  assumptions  usually  adopted  to  obtain  the 
set  of  primitive  equations  in  1-dimensions.  Then  the  method 
of  zonal-mean  averaging  is  outlined  to  show  how  the  set  of 
equations  is  reduced  to  2-dimensions. 

The  local  concentration  of  a  trace  gas  is  governed  by 
the  three-dimensional  continuity  equation 

(-£r  +  v  •  V)  f  =  Q/p  (7.2-1) 

d  t 

where  f(t,x)  is  the  mixing  ratio,  v(t,x)  the  velocity  wind 
fields  describing  the  general  circulation,  p(t,x)  the  air 
number  density  and  Q(t,x)  is  the  local  net  production  or 
loss  (by  chemical  and/or  physical  transformation)  of  the 
trace  gas.  Equation  (7.2-1)  gives  the  time  rata;  of  change 
of  f  in  the  Eulerian  description  of  fluid  motion.  The  quan¬ 
tities  f,  v,  p  and  Q  are  to  be  considered  as  Eulerian  field 
quantities  as  functions  of  time  and  spatial  location  wit  It 
coordinates  x. 
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In  order  to  solve  equation  (7.2-1)  for  the  species  con¬ 
centration,  one  must  be  able  to  provide  values  of  v  and  tem¬ 
perature  T  (T  is  necessary  for  calculation  of  reaction  rates) 
as  functions  of  space  and  time  either  by  parameterization  or 
by  solving  the  system  of  dynamic  equations.  The  atmospheric- 
circulation  is  governed  by  the  coupled  system  of  dynamic  and 
thermodynamic  equations  (cf.  Lorentz,  1967): 


momentum 

dy  = 

dt 


equation 
-2(2  xv- 


V4> 


thermodynamic  equation 


de 

dt 


J 


(7.2-2) 


(7.2-3) 


continuity  equation 

at  -  -p»-i  ,7-2-41 


equation  of  state:  ideal  gas  law 

p  =  p  R  T  (7.2-5) 

where  -3—  =  —  +  v  -V  is  the  total  time  derivative, 
dt  3 1 

The  above  are  to  be  considered  as  equations  for  the 
Eulerian  field  variables  v,  p,  p  and  0.  The  newly  intro¬ 
duced  symbols  have  the  following  meaning: 

S2  =  angular  velocity  of  earth 
M  =  average  mass  of  an  air  molecule 
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p 


pressure 


<>  -  the  geopotential  gz  where  g  is  acceleration  due 

to  gravity,  z  is  geometrical  altitude 

6  =  the  potential  temperature  related  to  temperature 

T  by  0  =  with  K  =  R/C  ;  k  the  gas  constant 

P  P 

and  C  the  specific  heat  at  constant  pressure 

T 

J  =  the  diabatic  influence  with  C  J—  the  heating 

p  u 

rate  per  unit  mass 


k  =  the  gas  constant 
T  -  temperature. 

Note  that  we  have  left  out  frictional  forces  in  the  momen¬ 
tum  equations  (7.2-2)  under  the  assumption  that  they  are 
unimportant  for  large  scale  motion. 

The  system  of  equations  (7.2-2)  to  (7.2-5)  is  coupled 
to  the  specie  equation  through  the  J  term  which  depends  on 
distribution  of  gases  such  as  ,  CO^,  N^O  and  Cll^  in  the 
atmosphere.  Thus,  in  principle,  equation  (7.2-1)  through 
(7.2-5)  must  be  solved  simultaneously  as  a  system. 

In  practice,  the  system  of  equations  presents  a  lormid- 
ubLe  numerical  problem  and  put  enormous  demand  on  both  com¬ 
puter  core  memory  and  computation  time.  This  is  particularly 
true  if  one  is  interested  in  a  realistic  chemical  scheme  in 
order  to  simulate  the  distribution  of  the  various  species 
in  the  atmosphere.  Besides,  the  set  of  exact  equations 
also  simulates  phenomena  of  little  interest  for  large  scale 
motions.  The  following  physical  assumptions  are  usually 
adopted : 
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A)  Replacement;  of  the  vertical  momentum  equation  by 
the  hydrostatic  equilibrium  condition,  i.e.,  the 
pressure  yradient  force  is  balanced  by  yeopotentral 
term) 


Tz  =  ~M  p  g  .  ( 7 . 2-> 6 ) 

it  is  observed  LhaL  motions  due  to  deviation  away 
from  hydrostatic  equilibrium  are  restricted  to 
oscillations  about  the  equilibrium  state  with 
time  scales  of  order  hour.  The  adoption  of  (7.2-G) 
effectively  filters  out  vertically  travelling 
sound  waves. 

B)  Discard  terms  containing  the  vertical  velocities 
in  the  remaining  two  components  of  the  momentum 
equation.  This  approximation  appears  to  be  justi¬ 
fied  because  of  the  fact  that  the  vertical  compon¬ 
ent  of  the  velocity  is  about  two  orders  of  magni¬ 
tude  smaller  than  the  horizontal  components. 

C)  Replace  r  in  the  resulting  equation  by  a  where  r 
is  the  distance  from  the  earth's  center;  a  is  the 
radius  of  the  earth. 

The  above  assumptions  lead  to  a  system  of  equations 
usually  referred  to  as  the  primitive  equations.  Next,  it 
would  be  desirable  to  write  the  momentum  equation  in  scaler 
form.  For  this  purpose,  it  is  convenient  to  introduce  a 


7-5 


7-6 


In  this  coordinate  system/  the  primitive  equations 


are 


In  this  arrangement,  equations  (7.2-10)  through  (7. 2-1.5) 
are  prognostic  equations  for  f,  u,  v  and  0  respectively, 
while  (7.2-14)  and  (7.2-19)  can  be  considered  as  diagnostic- 
equations  for  w  and  p.  Since  p  is  beiny  used  as  an  indepen¬ 
dent  variable,  equation  (7.2-10)  serves  the  purpose  of  trans¬ 
forming  between  geometrical  altitude  and  the  pressure  coor¬ 
dinate  . 

In  the  zonal  mean  model,  the  set  ol  equations  is  reduced 
to  two  dimensions  by  the  Eulerian  mean  averaging.  The  Euler- 
ian  mean  operator  averages  over  one  of  the  coordinates  while 
holding  all  other  coordinates  fixed.  In  the  case  of  the 
atmospheric  models,  the  integration  over  the  longitude  is 
done  along  constant  latitude  circles,  altitudes  and  time. 

Thus,  given  y  ( L  ,  •> ,  p ,  A )  ,  (where  t  is  time,  ;  the  latitude, 
p  the  pressure  height  coordinates  and  A  the  longitude) ,  one 
obtains 


y  (t ,  i|>,p) 


A  q  ( t  ,  y  ,  p  ,  a  )  d  A 

/  dA 


(7.2-17) 


It  is  convenient  to  define  g',  tne  deviation  from  the 
mean,  by 


g'  (L,<l,p,A)  -  g  (t ,  v  ,p)  -  y  ( t ,  v  ,  p ,  a  ) 


(7.2-18) 


i  l-  follows  from  (7.2-13  that 

y'  (t,(j>,p,A)  =  0 


(7.2-19) 
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and 


gh  =  gh  +  g'h* 


(7.2-1 9) 


When  the  Eulerian  mean  operation  is  applied  to  (7.2-10) 
to  (7.2-16)  using  the  fact  that  commutes  with  the 

coordinate  differential  operators,  we  have 


3f  ^  1  3  #7  -  . .  A  3  -. 

---  +  - —  T-r  f  V  COSi  +  V—  (  1  W ) 

dt  a  COS(|l  d  (J)  N  dp 


-  V" 


(7.2-20) 


3u  1  _  __ 

_ 2  3  0 


a  cos  4) 

-  2ft  v  cos (p  =  -F 


3  (u  v  cos^)  +  (uw) 
dp 


(7.2-21) 


-2  4 

u  tan4> 


3v  1  3-2  , .  ,  3  .- 

3t  +  aTcos^j  3 *  (v  cos4,)  +  3p  <v  v/)  +  a 

+  2  u  ft  sintf)  -  ^4t  =  -F 

a  d<|i  v 


(7.2-22) 


30  .  1  3  ,-  -  3 

_  +  - -  — —  0  VCOSi)  +  -r—  (0  w) 

3t  a  cos(p  3(f)  3p 


J  -  F. 


(7.2-23) 


1  3  .  3w 

— — — r  -rrr  (v  COS^/)  +  r—  =  0 

a  cos 4'  34>  dp 


(7.2-21) 


P  =  u  -J. 


(7. 2-25) 


3  z  _ _ 1_ 

9  3p  “  Mh 


R  T 
M  p 


(7. 2  —  2 1> ) 
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where  F 


- .  r  (  1  '  v  *  ee  s  v )  <  -  t 1  '  w  '  ) 

d  COS0  d P>  ap 


1 


-=-7-  (v'u'  COS  0)  +  -r—  (u'w1) 

2  ,  Dij)  T  up 


v 


0 


(v'^cos4>)  +  (v'w1)  i 


d  COS  if 

_ ]_ _ 0_ 

a  cos<f  3  0 

- - - 7  (  0  '  V  '  cost?)  +  7—  (0'w') 

a  cos(|>  3  0  op 


<) 

3p 


u ' 


l  an  v 


are  the  eddy  flux  divergence  terms.  Note  that  in  (7.2-21) 


u  v  tany 


)  into  the 


we  have  incorporated  the  curvature  term  ( 

term.  It  should  be  emphasized  that  the  eddy  flux  terms 
are  not  constrained  by  the  system  of  equations  and  must  be 
specified  or  parameterized  in  terms  of  other  variables. 

In  zonal-mean  models,  a  further  physical  assumption  is 
usually  adopted.  From  scale  analysis,  it  can  tie  shown  (cf. 
Iiolten,  1975)  that  equation  ( 7  .  2-2 2 )  reduces  to 


2Q  sin0  u  = 


3  JLs 

3  0 


(7. 2-27) 


the  lowest  order  of  a  Rossby  number  expansion.  This  is  the 
geostrophic  assumption  which  provides  a  good  approximation 
for  midlatitude  dynamics  though  some  authors  (Harwood  and 
Pyle,  1975)  argued  that  the  approximation  is  reasonable  up 
l-o  about  5°.  'Die  geostrophic  assumption  filters  out  occur¬ 
rence  of  gravity  waves  and  other  equatorial  waves  that  cause 
tin;  quasi-biennial  oscillations  in  the  stratosphere.  (1.2-21) 
is  awkward  to  use  as  z  is  not  one  ot  the  dependent  variables, 
one  can  differentiate  the  equation  with  respect  to  p  and  use 

the  hydrostatic  equation  to  obtain  the  thermal  wind  relation 
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(7.2-28) 


3u  _  -  _M_  Vpn  ’  1  a_o 

^  8p  2fi  sin<{>  a  3<J> 

to  replace  (7.2-22)  .  This  eliminates  Fv  as  an  extra  vari¬ 
able.  However,  one  still  needs  ways  of  determin ing  F f  ,  F 
and  FQ  to  close  the  system  of  equations. 

All  existing  zonal-mean  models  use  the  method  of  param¬ 
eterization  via  eddy  diffusion  tensor  originally  proposed  by 
Reed  and  German  (1965) .  Using  mixing  length  type  arguments, 
Reed  and  German  argued  that  fo^.  a  quantity  x  which  is  con¬ 
served  along  its  flow,  the  eddy  fluxes  are  related  to  the 
gradient  zonal  mean  quantity  x  via 


where  K  is  a  symmetric  tensor.  Following  the  procedure  sim¬ 
ilar  to  the  one  suggested  by  Reed  and  German,  Luther  (1973) 
analyzed  the  heat  transfer,  temperature  and  wind  variance 
data  of  Oort  and  Rasmussen  (1971)  and  derived  a  set  of  K 
tensors  as  a  function  of  time  and  space.  It  should  be  noted 
that  the  data  only  covered  part  of  the  northern  hemisphere, 
extrapolations  were  used  to  deduce  the  K's  at  places  where 

there  is  no  data  based  on  results  of  Newell  et  al.  (1966) 
and  Wofsy  end  McElroy  (1973).  Furthermore,  the  K's  for  the 
southern  hemisphere  are  obtained  by  reflecting  the  northern 
Hemispheric  values  in  the  appropriate  seasons.  The  set  of 
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r 


K  '  L  rout  LuLlior  (  1  *)  7  3 )  is  by  l.u  the  mu l  dimpleti-  wiiii  a 
set  of  values  for  each  month  covering  the  meridional  plane 
from  the  surface  to  Oil  km.  in  almost  all  oi  the  zonal 
models,  this  set  of  K1  s  provides  the  basis  l’oi  eddy  trans¬ 
port.  We  would  like  to  emphasize  here  Lhe  limrtaLion  oi 
this  approach: 

(i)  the  original  argument  of  Reed  and  German  (19G5) 
was  presented  to  treat  turbulence  diffusive  type 
motions.  Thus,  the  approach  may  not  be  appropri¬ 
ate  if  there  is  organized  wave-type  motion  in  the 
zonal  direction.  This  has  been  demonstrated  in 
studies  of  3-D  circulation  by  Mahlman  (1975), 
Matsuno  (1975) ,  Matsuno  and  Nakamura  (1979)  and 
Plumb  (19/9)  in  the  case  oi  planetary  wave  motions. 
To  date,  no  satisfactory  justification  for  apply¬ 
ing  the  eddy  diffusion  theory  to  stratospheric 
motion  has  been  presented. 

(ii)  The  argument  only  applies  to  conservative  flow. 
Thus,  strictly  speaking,  it  can  only  be  applied 
to  inert  tracers  or  to  potential  temperatures  in 
adiabatic  flow.  However,  in  zonal-mean  models, 

it  is  applied  to  chemical  species  with  finite  chem¬ 
ical  lifetimes  and  to  the  potential  temperature 
where  the  condition  of  adiabatic  flow  is  not 
strictly  satisfied. 

The  treatment  of  F  presents  special  problems  and  will  be 
discussed  later. 
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Next,  we  will  examine  how  to  handle  the  equations  lor 
determining  u,  v  and  w. 

Equations  (7.2-21)  and  (7.2-23)  can  be  written  as 


If  -A 


(7 . 2-30) 


=  (h 

at 


(7.2-31 ) 


where  a.nd$  do  not  contain  any  time  derivatives.  If  one 
differentiates  (7.2-30)  by  p  and  (7.2-31)  by  4>  and  takes 
the  appropriate  linear  combination  of  the  resulting  equa¬ 
tion,  by  virtue  of  the  thermal  wind  relations  (7.2-28)  one 
obtains 


d  ,  &<g.>K  130  ,  „ 

9p  2  \l  sin^  a  3  <p 


(7.2-32) 


Note  that  there  is  no  time  derivative  in  equation  (7.2-32). 
Equation  (7.2-24)  implies  the  existence  of  a  function  ip  where 


v  =  - 


w  - 


1  Dip 
cosep  9p 

1  .  it 

a  cos<p  9 cp 


(7.2-33) 


Using  (7.2-33)  to  illiminate  v  and  w  in  (7.2-32),  We  obtain 
a  second  order  partial  linear  differential  equation  for  ip 
where  the  coefficients  are  functions  of  u,  0,  J,  F^,  F  and 
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their  spatial  cicrj.vuLj.veii.  The  systeai  oi  equations  can  now 
be  solved  as  follows.  Given  f,  0,  u,  v  and  w  at  one  instance 
rn  time,  one  can  compute  Q/p,  J,  F^,  Fq  and  F^  and  solve  the 
prognostic  equations  (7.2-20),  (7.2-21)  and  (7.2-23)  for  f  >  u  an<3 
0  at  a  later  time.  With  these  new  values,  one  solves  (7.2-32) 
for  if  with  appropriate  boundary  conditins  and  generate  new 
v  and  w  via  equation  (7.2-33). 

We  will  next  discuss  how  to  obtain  Q/ p  once  a  photo¬ 
chemical  scheme  has  been  set  up.  The  term  Q/p  consists  of 
sums  of  terms  each  of  which  take  one  of  the  following  forms 


p  k.  . <T)  f . f  .  ; 
JO  id 


*2|‘iJk‘T>Vjfk 


for  photolytic  processes  and  two-body  and  three-body  reac¬ 
tions  respectively.  Note  that  the  temperature  dependence 

of  l.l»e  reaction  rates  k .  •  and  k.  ..  are  explicit  Ly  drspluyed. 

ij  J)K 

The  Longitudinal  behavior  o£  ^ ^  and  f^  arrse  from  the  diurnal 
effect  since  each  longitude  is  at  a  different  local  time  at 
any  instant  rn  Lime.  A  similar  problem  arises  in  the  1-b 
model  and  ls  usually  taken  Care  of  by  putting  diurnal  aver¬ 
aged  photolysis  rates  in  the  chemical  scheme  so  that  the 
resulting  species  concentrations  represent  diurnally  averaged 
values.  in  addition,  the  longitudinal  distribution  of  the 
gas  may  be  affected  by  wave  motion  in  the  zonal  direction. 

No  attempts  are  made  to  treat  such  effects.  in  the  zonal 
averaging  process,  the  term  Q/p  is  obtained  by  taking  the 
sum  of  the  corresponding  products  of  the  zonal-meaned  quan¬ 
tities.  Eddy  flux  terms  are  ignored.  The  treatment  of  d  is 
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similar  in  that  the  term  is  obtained  by  using  zonal  ly-iueaned 
values  in  the  expression  while  ignoring  eddy  flux  terms. 

The  treatment  of  i-'u  is  more  diiliculL  because  the  u 
momentum  is  definitely  not  conserved  in  the  flow.  Attempts 
to  parameterize  in  terms  of  potential  vorticity  and  poten¬ 
tial  temperature  have  had  little  success  (Green,  1970;  Wiin- 
Neilsen  and  Sela,  1971).  Vupputuri  (1979)  employed  an  ad 
hoc  parameterization  while  Harwood  and  Pyle  (1975)  deduced 
the  momentum  fluxes  from  Nimbus  V  SCR  data.  Until  a  more 
satisfactory  theoretical  approach  is  available,  the  satel¬ 
lite  data  approach  seems  to  be  preferrable  at  this  stage. 

7.3  AER  2-D  Model:  Intermediate  Results 

During  the  second  stage  of  the  development,  we  have 
concentrated  our  efforts  on  interfacing  the  chemical  scheme 
with  the  rest  of  the  model.  At  the  present  stage,  the  model 
only  solves  (7.2-20)  for  the  diffusive  species.  The  values 
necessary  for  v,  w,  K's  are  read  in  from  input  files.  A 
description  of  the  model  can  be  found  in  Appendix  B. 

Three  series  of  chemical  studies  were  made  using  the 
oxygen  chemistry,  the  oxygen  +  HOx  chemistry  and  the  full 
O  -HO  -NO  -Cl  chemistry.  In  each  of  the  studies,  we  have 
concentrated  on  the  latitudinal  and  seasonal  behavior  of  the 
calculated  ozone  column  densities.  This  approach  is  taken 
because  0^  is  a  dominant  factor  in  determining  the  behavior 
of  other  trace  gases.  In  addition,  observations  of  the  ozone 
network  over  the  past  10  years  have  provided  a  data  base  will) 
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complete*  global  coverage  appropriate  for  model  validation 
analysis.  Fig.  7.1  shows  the  observed  ozone  column  density 
as  a  function  of  latitude  and  season.  The  calculated  results 
will  be  compared  with  Fig.  7.1  to  judge  how  well  they  agree 
with  observations . 

In  the  oxygen  chemistry,  the  classical  Chapman  scheme 
is  adopted 


C>2  +  hv  -*■  0  +  0  J2 

03  +  hv  ■*  O2  +  0  J3 

02  +  0  +  M  -<•  03  +  M  k2 

O  +  03  -*■  202  k3 


(7.3-1) 


where  the  J's  are  photolysis  rates,  k's  are  reaction  rate 
constants . 

The  production  term  takes  the  form 


Qi  2j2l02]  ”  k2[02] [M] 


(7.3-2) 


Taking  [0^1  =  [Ml/4.75,  equation  (7.2-20)  is  solved  with  the 
above  production  terms.  The  calculated  C>3  column  densities 
as  a  function  of  latitude  and  time  are  given  in  Figure  7.2. 
Note  that  the  overall  features  of  the  spatial  and  temporal 
behavior  of  the  observed  G3  distribution  are  reproduced. 
These  include: 

The  spring  time  column  ozone  maximum  in  the  northern 
hemisphere  at  70°  N  and  extending  all  the  way  to  the 
pole  . 
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OBSERVED  OZONE  COLUMN 
(Dobson  Unit) 


Figure  7 . 1 
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MONTH 

OZONE  COLUMN  (Dobson  Unit) 
(Ox) 


Figure*  1.2 

Calculated  ozone  column  densely  as  a  Lunction 
ol  latitude  and  season  using  oxygen  chemistry  only 


.  The  spring  time  column  ozone  maximum  in  Lhe  southern 
hemisphere  at  about  60°  S. 

.  The  relatively  small  seasonal  variation  ol  the  ozone 
column  in  the  equatorial  region. 

.  The  relative  magnitudes  of  the  two  maxima  at  high 
latitudes  compared  with  the  values  at  the  equatorial 
region. 

However,  the  magnitude  of  the  column  densities  is  about  a 
factor  of  2  larger  than  the  observed  values,  suggesting  that 
certain  loss  terms  had  been  ignored. 

In  the  next  series  of  studies,  the  O  +  HO  chemistry 
is  employed.  The  reactions  included  are 


H20  +  0(1D) 

2  OH 

ka 
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-* 
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OH 
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The  chemical  production  term  for  takes  the  form 


7-19 


2 J2 102] 


(7.3 


Qi  = 


-  2k3  101  1 0 3 J 


kf[ou][o3]  -  kinio2i[o3l  -  2k  [O]  I  no. 


Thu  calculated  ozone  column  densities  are  shown  in  Fig.  7.3. 
The  magnitude  of  the  column  densities  is  reduced  by  a  factor 
of  1.8  compared  to  the  oxygen  case.  it  is  interesting  to 
note,  however,  the  latitudinal  and  seasonal  behavior  remains 
remarkably  similar  to  the  oxygen  case.  This  seems  to  sug¬ 
gest  that  although  local  photochemistry  may  be  responsible 
for  the  magnitudes  of  the  calculated  ozone  column,  dynamical 
transport  plays  an  important  role  in  determining  the  corres¬ 
ponding  seasonal  and  latitudinal  behavior. 

in  simulating  the  present  atmosphere,  the  time  depen¬ 
dent  diffusion  aquations  are  solved  for  the  following  species 
Cix  und  its  precursor  molecules  via  Cll-jCl ,  CU3CC13. 
CC14,  F- 11,  F- 1 2 

N0X  and  its  precursor  molecule  N2o 

.  o3 . 

Ln  addition,  the  following  species  are  calculated  assuming 
photochemica 1  equilibrium: 

.  0  (XD)  ,  0( JP) 

.  H,  OH,  H02,  H202 
.  Cl,  CIO,  HC1,  H0C1 ,  ClNOj 
NO,  NO.,,  UNO  3 . 

fixed  mixing  ratio  boundary  conditions  are  used  on  all 
the  diffusion  species  except  P-1 l  and  F-12  for  which  time- 
dependent  flux  boundary  conditions  are  used.  The  release  of 
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LATITUDE 


OZONE  COLUMN  (Dobson  Unit) 
(Ox  +  HO* ) 


Figure  7.3 

Calculated  ozone  column  densities  as  a  function  of 
latitude  and  season  using  {oxygen  +  HO  )  chemistry. 

A 


l’-ll  and  F-12  is  restricted  to  between  2  5"  N  to  60°  N.  The 
release  rate  follows  the  data  compiled  by  the  CMA  luoro 
carbon  Panel  (1980).  in  the  present  model  calculation,  CU^, 

CO  and  H^O  are  specified  as  a  function  of  latitude,  altitude 
and  t ime . 

The  distributions  of  NO^  and  Cl^  play  an  important  role 
in  controlling  the  local  radical  concentrations  which  in  turn 
determine  the  local  ozone  production  and  loss  rates.  figure 
7.4  shows  a  typical  calculated  distribution  of  Clx  from  our 
simulation.  It  is  worthy  of  note  that  in  the  case  of  C 1 x  and 
NO,.,  the  concentrations  at  the  equator  are  lower  than  the 
corresponding  concentrations  at  high  latitudes  m  the  lower 
stratosphere.  This  is  indicative  of  the  upweliing  air  motion 
at  the  equator  and  subsequent  transport  of  material  towards 
higher  latitude.  In  our  calculation,  the  upper  stratosphere 
concentrations  of  Cl  and  NO  are  2  ppbv  and  17  ppbv  respec- 

X  X 

t  .1 '/  e  L  y  . 

figure  7 .  ‘j  shows  the  latitudinal  and  seasonal  variation 
or  calculated  ozone  column.  The  magnitude  of  the  ozone  col¬ 
umn  compares  favorably  with  observation  around  the  equatorial 
region.  However,  the  magnitude  of  the  spring  time  maxima  in 
both  hemispheres  is  about  30%  too  low  compared  to  observation 

Figure  7.6  shows  the  number  density  distribution  of 
ozone  corresponding  to  the  month  of  April.  The  maximum  num¬ 
ber  density  occurs  at  25  km.  The  slightly  higher  concentra¬ 
tion  near  the  North  Pole  is  indicative  of  the  spring  time 
maximum.  Again,  compared  to  observation  (cf.  Outsch,  1978), 
i  in-  calculated  mini  bet  density  around  the  North  i’ole  is  on  the 
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LATITUDE 


ax  MIXING  RATIO  (ppbv) 
APRIL 

Figure  7.4 

Calculated  present-day  Clx  mixing  ratio  as  a  function 
of  altitude  and  latitude  for  the  month  of  April. 
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LATITUDE 


OZONE  NUMBER  DENSITY  (nb) 

April 

Figure  7.6 

Calculated  present-day  ozone  number  density  as  a  function 
of  altitude  and  latitude  for  the  month  of  April. 


low  side.  Wo  will  be  performing  analyses  to  try  to  deter¬ 
mine  if  the  discrepancy  could  be  attributed  to  chemistry  or 
dynamic  transport. 

We  also  performed  a  sensitivity  study  on  the  response 
of  the  model  to  changes  in  Clx  concentration.  More  specif¬ 
ically,  we  solved  the  0^  diffusion  equation  using  twice  the 
Clx  concentration  of  the  present  atmosphere,  in  effect  adding 
an  extra  2  ppb  of  Clx  to  the  stratosphere.  The  result  shows 
that  this  leads  to  a  7%  reduction  of  the  ozone  content  of  the 
atmosphere.  The  %  column  reduction  as  a  function  of  season 
and  latitude  is  given  in  Figure  7.7.  Note  that  the  largest 
%  reduction  occurs  in  the  winter  hemisphere  at  high  latitudes 
while  the  smallest  reduction  occurs  around  the  equator.  The 
results  are  given  in  different  format  in  Figures  7.8a  and 
7.8b  which  clearly  show  the  seasonal  variation  at  fixed  lat¬ 
itude.  These  features  are  in  qualitative  agreement  with 
those  of  Pyle  (1978)  and  Vupputuri  (1979) . 

Figures  7.9  through  7.11  shows  the  altitude  and  latitude 
cross-section  of  %  ozone  reduction  for  the  months  of  January, 
April  and  July.  Figures  7.9  and  7.11  again  illustrate  that 
the  maximum  ozone  reduction  occurs  in  the  winter  hemisphere. 
Also  evident  from  the  figures  is  the  slight  local  ozone  in¬ 
crease  at  the  equator  around  10  to  20  km.  This  is  usually 
attributed  to  the  so-called  "self-healing"  effect  due  to  the 
increase  of  solar  flux  as  a  result  of  reduction  above  20  km. 
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LATITUDE 


OZONE  COLUMN 

%  OZONE 

(Dobson  Units) 

COLUMN  REDUCTION 

,  J  ,  F  ,  M  ,  A  ,  M  , J,JiAiSiOiNiPi 

MONTH 

Figure  7.8a 

Calculated  present-day  ozone  column  densities  and  %  ozone 
reduction  with  Clx  doubling  as  a  function  of  season  for 
fixed  latitude  in  the  "crthern  hemisphere. 


OZONE  COLUMN 
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%  OZONE  REDUCTION 


JANUARY 


Figure  7.9 

Calculated  %  ozone  reduction  for  Clx  doubling 
as  a  function  of  altitude  and  latitude 
for  the  month  of  January. 
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Figure  7.10 


Calculated  %  ozone  reduction  for  Clx  doubling 
as  a  function  of  altitude  and  latitude 
for  the  month  of  April. 
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Figure  7.11 

Calculated  %  ozone  reduction  for  Clx  doubling 
as  a  function  of  altitude  and  latitude 
for  the  month  of  July. 
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DOCUMENTATION  OF  THE  AER  I-D 
PHOTOCHEMICAL  DIFFUSION  MODEL 
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II.  Chemical  Scheme 
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IV.  Calculation  of  Production  and  Loss  Term 

V.  Treatment  of  Photochemical  Species 

VI.  Treatment  of  Diffusive  Species 

VII.  OUTPUT 

VIIL  Running  Instructions 
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I .  GENERAL  STRUCTURE 

The  AER  one-dimensional  photochemical  diffusive  model 
calculates  the  number  densities  of  chemical  species  in  the 
earth's  atmosphere  with  a  variety  of  user-selected  options. 

The  basic  design  of  the  program  has  been  discussed  in 
Section  3.  in  Figure  A-l,  a  more  detailed  flow  chart  of  the 
model  is  presented.  In  what  follows,  each  component  of  the 
model  will  be  briefly  discussed. 

In  this  section,  some  of  the  variables  used  in  the  pro¬ 
gram  will  be  defined.  All  variables  that  are  essential  to 
input,  output  are  included  in  the  following  discussion. 

The  atmosphere  is  divided  into  81  levels  each  1  km  thick. 
The  bottom  layer  is  at  ground  level.  All  altitude-dependent 
variables  are  given  in  terms  of  81  element  arrays;  these  in¬ 
clude: 

TIALT (81) :  Altitude  array,  gives  the  corresponding 

geometrical  altitude  of  each  level  in  cm. 

TITK(81):  Gives  the  temperature  of  each  level  (°K) 

EDC(81):  Gives  the  values  of  the  eddy  diffusion  co- 

2 

efficient  at  each  level  (cm  /sec) 

DM  (81):  The  air  density  at  each  level  (— --ec^*e.s, 

cm3 

SS  (3240)  :  Array  for  species  number  densities  (5!2ieculesj 

cm 3 

Number  densities  are  stored  in  40  records, 
each  81  elements  long.  The  density  of  the 
I  species  at  the  level  is  given  by 

SS( ( I- I)  X  81  +  J) 
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SM ( 1458 ) : 


Array  for  storing  tho  calculated  diurnal- 
averaged  member  densities  (n  (t , z ) )  for 
photochemical  species.  Number  densities 
are  stored  in  18  records  each  81  elements 
long.  The  mean  density  for  the  iL*‘  specie 
at  the  J*"*1  level  is  given  by  SM((1-1)  x  81  +  J) 
SPAIR(9f71)  :  Array  for  storing  the  calculated  diurnal- 
averaged  product  of  number  densities 
(n^(t,z)  xn^(t,z)).  These  terms  are  nec¬ 
essary  for  computing  the  production  and 
loss  rates  of  the  diffusion  species  in  the 
diurnal  mode 

SJ(2835):  Array  for  storing  calculated  photolysis 

rates  as  functions  of  altitude.  There  are 
35  records  each  81  elements  long.  The  pho¬ 
tolysis  rate  for  the  Ith  photolytic  reac¬ 
tion  at  each  level  J  is  given  by 
SJ( (1-1)  x  81  +  J) 

SR (4050):  Array  for  storing  the  calculated  altitude- 

dependent  reaction  rates.  There  are  50 
records  each  81  elements  long.  The  Ith 
altitude-dependent  reaction  rate  at  the  Jth 
level  is  given  by  SR((l-l)  x  81  +  j) 

Since  the  photochemical  species  are  solved  one  altitude 
at  a  time,  it  is  convenient  to  define  working  arrays  for 
variables  at  one  altitude;  these  include: 


SBUF(378):  Working  array  for  storing  the  time-dependent 

number  densities  at  the  altitude.  There  are 
18  records  each  21  elements  long.  The  num¬ 
ber  densities  for  the  It^>  specie  at  the 
time  step  is  given  by  SBUF((I-1)  x 21  +  K) 
[use  in  diurnal  option  only] 

J(50):  Working  array  for  storing  photolysis  rates 

at  the  altitude 

DSJ (580) :  Working  array  for  storing  the  time-dependent 

photolysis  rate  at  the  altitude.  There  are 
27  records  each  21  elements  long.  The  Ith 
photolytic  rate  at  Kth  time  step  is  given 
by  DSJ ( ( I- 1 )  x  21  +  K) 

R(75) :  Working  array  for  storing  all  reaction  rates 

at  the  altitude 

The  following  arrays  are  used  in  photolysis  calculation. 

The  array  index  corresponds  to  the  wavelength  dependence: 

o 

WL (72):  Gives  the  wavelength  of  the  index  in  A 

FLX(72):  Gives  the  unattenuated  solar  flux  integrated 

over  the  corresponding  wavelength  window 

( — ^ - ) 

cm*  sec 

QA(252Q) :  Gives  the  absorption  cross-section  (cm^)  for 

the  photolytic  processes.  There  are  35 
records  each  72  elements  long.  The  cross- 
section  of  the  Ith  process  at  the  Jth  wave¬ 
length  is  given  by  QS((I-1)  x 72  +  J) 

In  addition  there  are: 


SID (50) : 
QID (25) : 
RR(4 ,75) 


BCB  (50)  : 


BCT (50 )  : 


DCLT (21) 


FSF (50) 
ICF (919) 


Alphanumeric  symbols  for  chemical  species 
Alphanumeric  symbols  for  photolytic  process 
Parameters  for  computing  reaction  rate  con¬ 
stants.  Each  reaction  is  parameterized  by 
4  variables. 

Bottom  boundary  values  for  the  species.  Twc 
formats  accepted.  IF  0  <  BCB(I)  <  1,  BCB(I) 
is  assumed  to  be  the  bottom  volume  mixing 
ratio.  Otherwise,  BCB (I)  is  assumed  to  be 


the  boundary  flux  in 


molecules 
cm 2  sec 


Top  boundary  condition.  The  comments  on 
formats  in  BCB  apply  here  also 
Gives  the  time  intervals  between  time  steps 
in  seconds 

Array  necessary  for  closure  families 


II.  CHEMICAL  SCHEME 


The  species  included  arc: 

Photochemical  species 

Oxygen  family:  O(^D),  O 

Hydrogen  family:  H,  OH,  H02,  H202,  CH20 
Chlorine  family:  Cl,  CIO,  HCi ,  ilOCl,  C1NC3 
Nitrogen  family:  NO,  N02,  HNO3,  NO3,  N205,  H02N02 
Diffusion  species 

n2o,  ch4,  h2,  CH3CI ,  cci4,  CFCI3 ,  cf2ci2,  cf2o,  CO, 

NOX,  C1X,  FX,  O, 


fixed  species 

h2o,  o2,  M 

The  species  NOX,  C1X  and  FX  correspond  to 

NOX  (odd  nitrogen  species)  =  NO  +  N02  +  HN03  +  N03  + 

2xn2o5  +  ho2no2  +  cino3 

C1X  (odd  chlorine  species)  =  Cl  +  CIO  +  HC1  +  H0C1  + 

cino3 

FX  (odd  fluorine  speices) . 

The  chemical  species  react  with  each  other  via  photol¬ 
ysis,  two-body  and  three-body  reactions.  The  reactions  in¬ 
cluded  in  the  model  are  listed  in  Table  A-l  while  the  photo- 
lytic  processes  are  given  in  Table  A-2.  The  meaning  of  the 
reaction  rate  parameters  is  discussed  in  IV.  The  reaction 
rates  quoted  here  follow  the  NASA  (1979)  recommendations. 

In  addition  to  chemical  losses,  some  species  are  re¬ 
moved  by  heterogeneous  removal  processes  of  rainout  and 
washout.  Heterogeneous  processes  are  included  in  the  loss 
terms  of  HN03,  HC1  and  H202«  The  process  is  modeled  by  a 
pseudo-first  order  rate  in  the  troposphere  with  an  effective 
lifetime  of  approximately  7  days  (cf.  Wofsy,  1976). 

Finally,  boundary  values  are  prescribed  for  the  diffu¬ 
sive  species  to  simulate  exchange  with  the  ground.  These 
.ire  given  in  Table  A- 3 . 


TABLE  A-l 


REACTION  RATE  PARAMETERS 


T  Y‘->C 

D&(V  1 

P4R2 

PARI 

1 

2. 

]  .AAF-1 3 

1 .07c *01 

0. 

3M.CH1CL=XX*H30»3 

p 

1  . 

3. is^- i n 

0. 

o . 

hpo ♦ o i n-3»9M , i 

1 

1  . 

1  .40P-1 n 

0. 

0. 

C  S4 ♦ 0 1 Q=CH20*0H«  1 

U 

1  . 

1  .oof-i  n 

0. 

r\ 

• 

H?*Oiri  =  OH.H,  1 

•! 

3. 

1  .*ftF-32 

-3.40F*02 

0. 

H*03*w=hO 3*m,2 

4 

3. 

l  .aop-i  n 

4.70F *02 

o. 

-i*oi=nH  +  03 . 3 

7 

3. 

1 .60?-! 3 

9.40F  *03 

0. 

3H*03=H03*0?.2 

a 

3. 

1.0  (-11) 

7.5G27 

0  . 

1h*H?O2=H20*h02.2 

3 

3. 

1  .40c-14 

S. HOF *03 

0. 

MO3«'0.1  =  0H*3o0?»2 

i  n 

P  . 

1  .nn^'-in 

3. 50F  *03 

0  . 

H03*0=Om«03. 3 

1  l 

1  . 

*.  fin- -i  3 

0. 

0. 

H02*N0=0H*M92. 1 

1  ? 

1  . 

4 .  nor- 1 ? 

0. 

0. 

M02  ♦  Cl  0  =  HOCL  *02. 1 

1  7 

1  . 

3.S0C--1  3 

0. 

0. 

H0-2*riO3  =  H?C13*n2. 1 

)4 

1  . 

4 » 0  0  F  - )  1 

0. 

0. 

h03*0m  =  h20*O2. 1 

i  s 

1  . 

4.S0F-1 1 

0. 

0. 

-I02*CL  =  HCL*0?»  1 

i  s 

1  . 

1  .OOF-1 1 

0. 

0. 

H ♦ HO 2= HP* 03* 1 

1  7 

3  . 

4.20F-1  o 

9  .  SOF  *  02 

0. 

H*H02=2u0h, p 

i  ^ 

3. 

3 .  i  o  f  - 1 1 

1 . ISP *03 

0. 

1h*C0  =  H*  XX  ,  3 

1  3 

3. 

1  .102-1 n 

3 . SOP  +  03 

0. 

1H*0  =  H*O2  .  3 

?r> 

3. 

3 . 0  0  r  -  1  3 

'3.  1  3  F  ♦  0  3 

0. 

DH.HCl  =H20*rL.? 

31 

3. 

2. 3HP-1 2 

1 ,71F*03 

0. 

1H*CH4=CH?0*XX  .2 

33 

1. 

?.«0 f*nn 

3 . 4  0  F  -  l  1 

1 . 30F ♦ on 

1H*t\J03  =  HN0l*M,4 

?i 

1  . 

s.snr-i a 

0. 

0. 

lw*HM03=H2n*N03. 1 

24 

1  . 

1  .70P-1 3 

0. 

0. 

1H*0H=H20*1. 1 

31 

3. 

*  .  7  0  ?  -  )  3 

3 . 0 1 F ♦ 0  7 

0  . 

Dh*H2=H30*H. ? 

?4 

3. 

S  .  /  0  'T  -  1  1 

2 . 40F  +  03 

0  . 

CL*H2=HCL*H. 2 

2  7 

1  . 

<S.nn=--i  i 

0. 

n . 

OL  +H203  =  HCL  *H02 • 1 

?« 

0  . 

1. 4  0  F  ♦  0  0 

1 .SOF-1 1 

1  .  9  0  E  ♦  0  0 

CL0*NO2=CL^O3*M.4 

33 

3. 

3.S0P-1  1 

2 •  57F  +  03 

0. 

01*03  =  0.0*02,2 

in 

3. 

p.'rof-i  3 

1 ,1SF*m 

0. 

CL  *CH4  =  hCL  *CH30. 3 

1 1 

3. 

7. 70^-1  1 

1  . 1 0  F  ♦  0  3 

0  . 

on*o=ci  *02.2 

i? 

1  . 

3 . 0  0  c  -  1  1 

0. 

0. 

0L0*NIO  =  CL*MO2.  I 

1 1 

1  . 

0. 

0 . 

0. 

0 ♦ wCL  =  OH ♦ Cl  .  1 

14 

3. 

1 .  osf-14 

-5. 30 F *02 

0. 

3*0?*w=03*^,2 

3  1 

p  . 

l .  ior-n 

3 . 1  4F  *01 

0. 

0*03  =  2^02, 3 

14 

1  . 

1  .ftOE-1 3 

0. 

0. 

M03*N03  =  M30Sf  1 

1  7 

3. 

3.1nr -) 2 

1 . 4SF  ♦  0 .1 

0. 

'10*03  =  ^02*03,2 

IQ 

3. 

1  . 30^-1 1 

3  •  45F  *01 

0. 

\|02*01  =  M03*02.3 

11 

1  . 

9.  1  2P-1  2 

0 . 

0  . 

vjO2*0  =  MO*03.  1 

A  n 

1  . 

5.  noE- ) 1 

0. 

0  . 

V20*  Dl 0=XX  *02*  1 

4  1 

1  . 

S  .  0  0  2  -  1  1 

0. 

0. 

m?o*oi  n=2«\io.  1 

U  P 

1  . 

1.00r-l  (1 

0  . 

0  . 

11  0*CPCI.3=xy  *CLO. 

4  ■> 

1  . 

l.nor-1 n 

0. 

o. 

1)  0*CP3ru3=xX*CL<i 

TABLE  A- 1 

(cont . ) 

Jj 

T  Y  ' J  P 

w  AW  1 

DAP? 

PARI 

t-i 

uu 

1  . 

o. 

0. 

0. 

91U«-Ch3CL  =  0hyXX,1  ’  '  ;  ! 

4  5 

)  . 

n. 

0. 

0. 

DlD»CrL4=XX»CL0. \ 

4«, 

?  . 

?, 1  OP-1  1 

-1 . 0  7P ♦ 0  ? 

0. 

010*8=0 *8,? 

4  / 

0  . 

5.00E+0O 

6.50E-1? 

5. 00f  *00 

-l02*\!O2  =  H0?\'02*M,4  r! 

ua 

1  . 

3 . 0  0  c  -  1  0 

0. 

0. 

OH  +  HOn. =h?o*clo. i 

1  . 

1 

0. 

0. 

01D*CP2O=XX . 1  , 

5 n 

1  . 

1  .  0  0  P  -  1  A 

0. 

0 . 

3  ♦  HOCl  =OH*  CL  0  »  1  '  1 

51 

1  . 

1  JnoE-14 

0. 

CL*HOCL=HCL*CLO, 1  U 

52 

1  . 

cv 

1 

u.' 

o 

• 

a 

0. 

0. 

CLO*  Oh=hO?*CL«  1  hj 

53 

2  ■ 

3.50F-1 P 

1  . SAP  *03 

0. 

C^oCCl  3*0h=yX*H?O,?  .j 

54 

1  . 

4  .  o  n  E  - 1 5 

0. 

0. 

021 03  =  2»0?*0. 1  1 

55 

p. 

0. 

0. 

0. 

CL0*HO2=HCl*03»? 

56 

p. 

0. 

0. 

0. 

CL0*0w=HCL*O2.2  h 

57 

i  . 

1  .  0  ()  p  - 1  1 

0  . 

0. 

ch?o+oh=5?o*co*h. i  : 

56 

?. 

9 . ?QP- ] 1 

6 • 90F  +  0  1 

n. 

CL*CH?0  =  HCl*CO*H,;>  V 

59 

p. 

5. OOE-1  ? 

2.00F+0? 

o. 

-•02NC>P*0H  =  HP0*NO?,p 

60 

l  . 

1 .00P-1 P 

0. 

0. 

HO?NDp*CL=HCL*NOP. 1  ; 

61 

p. 

5. P0r-0* 

1  .01  F-0'f 

o. 

402NDP*m=hOP*N02.? 

62 

p. 

1 .OOF-15 

0. 

0. 

Hf>2NO?*r.LO  =  HOCL*NO?*  p  |] 

63 

?. 

3. OOF-  1  ? 

8.08F+02 

0  . 

H0?N0?*O=0H*NO2.P  jj 

64 

p. 

3. OOF-1 ? 

8.08F*0? 

o. 

CLki03*0  =  CL9*ND3.P  li 

TABLE  A- 2 


PHOTOLYTIC  PROCESSES 


°2 

+ 

hv 

-** 

O  +  O 

°3 

+ 

hv 

-> 

o2  +  0 

°3 

+ 

hv 

-b 

o2  +  0(LD) 

no2 

+ 

hv 

-*■ 

NO  +  0 

hno3 

+ 

hv 

-b 

OH  +  N02 

no3 

+ 

hv 

-b 

NO  +  02 

N2°5 

+ 

hv 

-b 

no2  +  no3 

n2° 

+ 

hv 

-b 

n2  +  o 

C1N03 

+ 

hv 

-b 

CIO  +  N02 

HC1 

+ 

hv 

-b 

H  +  Cl 

HOC1 

+ 

hv 

-> 

OH  +  Cl 

^2^2 

+ 

hv 

-► 

OH  +  OH 

ch2o 

+ 

hv 

-► 

2H  +  CO 

ch2o 

+ 

hv 

h2  +  CO 

h2° 

+ 

hv 

-b 

H  +  OH 

ho2no2 

+ 

hv 

-b 

HO 2  +  N02 

cfci3 

+ 

hv 

-b 

product 

CF2C12 

+ 

hv 

->■ 

product 

CC14 

+ 

hv 

-> 

product 

CH3C1 

+ 

hv 

-► 

products 

cf20 

+ 

hv 

-+• 

products 

CH3c13 

+ 

hv 

-b 

products 

It  is  assumed  that  CFCl^,  CF2C12,  CCl^,  CH-jCl  and  CH^Cl^ 


release  all  their  chlorine  atoms  following  initial  reactions. 


TABLE  A- 3 


BOUNDARY  CONDITIONS 


Fixed  Mixing  Ratio  Boundary  Conditions 


n2o 

300 

ppbv 

CH  . 

4 

1.5 

ppmv 

H2 

500 

ppbv 

ch3ci 

1.2 

ppbv 

cci4 

130 

pptv 

CO 

200 

ppbv 

NOX 

100 

pptv 

C1X 

1 

ppbv 

°3 

18 

ppbv 

Flux 

Boundary  Conditions 

CFCI3 

9.3xl05 

molecules 
cm2  sec 

cf2ci2 

1. 3xl06 

molecules 
cm2  sec 

ch3cci3 

l.OxlO7 

molecules 
cm2  sec 

cf2o 

1 . 0x1 03 

molecules 
cm2  sec 

FX 

1. 0xl03 

molecules 

cm2  sec 
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Ill . 


INPUT 


The  BLOCK  DATA  contains  the  information  necessary  to 
initialize  the  variables  in  the  program..  Any  changes  can 
be  made  directly  in  the  routine.  The  block  data  values  are 
overridden  by  the  use  of  Namelist  INPUT,  or  tape  input. 

BLOCK  DATA: 

The  BLOCK  DATA  is  separated  into  4  routines  BDl ,  BD2 , 

BD3  and  BD4 . 

BDl:  Contains  a  set  of  values  for  all  chemical  species 

concentrations.  These  correspond  to  the  values 
used  for  the  fixed  species  throughout  the  program 
and  serve  as  first  guess  values  for  the  calculated 
species,  BDl  also  contains  an  identity  array  for 
all  the  chemical  species  giving  the  chemical  sym¬ 
bol  of  the  specie  in  alphanumeric  format.  This 
facilitates  an  easy  search  procedure.  The  BLOCK 
DATA  values  can  be  overridden  if  one  decides  to 
read  in  previously  calculated  number  densities 
attached  as  an  input  file  (tape  7  input)  . 

BD2:  Contains  the  necessary  information  for  setting  up 

the  chemical  reaction  scheme.  It  consists  of  an 
array,  RR{4,75),  specifying  4  parameters  for  cal¬ 
culating  the  reaction  rate  constant  for  each  of 
the  reactions.  It  also  contains  the  absorption 
cross-sections  for  photolytic  processes.  Both 
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RR  and  QS  can  be  entered  via  namelist  to  override 
BLOCK  DATA  values. 

BD3:  Contains  the  standard  atmospheric  temperature  as 

a  function  of  altitude;  the  array  giving  the  inte¬ 
grated  unattenuated  solar  flux  as  a  function  of 
wavelength  intervals. 

BD4 ;  Contains  parameters  necessary  for  computing  the 
eddy  diffusion  coefficients. 


ZB  =  altitude  of  bottom  boundary  (cm) 

ZT  =  altitude  of  top  boundary  (cm) 

H(cm)  =  the  thickness  of  each  interval;  program 
set  up  for  H  =  1.E5  cm  or  2.E5  cm 
MODE  =  0;  use  BLOCK  DATA  values  for  species  den¬ 

sities 


=  1;  read  species  density  profiels  from  a 

local  file  tape  7,  which  contains  pre¬ 
viously  calculated  number  densities 


default  value  is  0. 


Parameters  Cor  Chemical  Scheme: 


RR 


QS 


PJ03 


WO 


array  containing  parameters  for  calcu¬ 
lating  reaction  rates  (default  to  BLOCK 
DATA  values) 

array  containing  absorption  cross-sec¬ 
tions  (default  to  BLOCK  DATA  values) 
parameter  for  photolysis  rate  of 
Oj  +  hv  -*■  02  +  0  for  A  >  4025  A 
sedimentation  velocity  (set  =  0  for  mol¬ 
ecules) 


FSF, 

ICF 


LNOX 


=  indices  for  closure  family 


G22,  factor  for  adjusting  the  reaction  rate 

G28 , 

G47  in  diurnal-average  mode 

GAMMA  =  factor  parameterizing  the  production  of 


O^  due  to  methane  oxidation 
=  altitude  index  for  NOX  injection 


ADDNOX  =  amount  of  NOX  injected  (m-0^cule^) 

J  cm3  sec 

TA  =  latitude  to  be  used  for  photolysis  calcu¬ 
lation  (default  =  30°) 

DA  =  declination  to  be  used  for  photolysis 
calculation  (default  =  0°) 

Parameter  for  Diurnal  vs  Diurnal-Average  Options 

SST  =  (T, F)  (diurnal-average  mode,  diurnal 


mode) 

NDIV  =  number  of  divisions  into  which  the  day¬ 
light  hours  are  divided  [Note,  some 
array  sizes  will  have  to  be  changed  if 
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ND IV  '•  6;  use  only  in  diurnal  model 
Boundary  Conditions: 

BCB  -  bottom  boundary  values 

BCT  =  top  boundary  values 

Convergence  Tolerance  Limit: 

EP1  =  convergence  criteria  for  photochemical 

species  lused  in  NEWTON) 

EP2  =  convergence  criteria  for  diffusive  species 
ITC  =  iteration  limits;  program  aborts  if  con¬ 

vergence  criteria  not  satisfied  after  ITC 
iterations 
PRINT  OUT  OPTION 

LEVP  =  level  of  print  out;  6  is  normal 
DZ  =  print  out  densities  at  intervals  of  DZ 

(cm) ;  default  =  2xl05  cm 
ZBP  =  print  out  density  profiles  from  height 
ZBP  (cm) ;  default  =  0 

ZTP  =  print  out  density  profiles  to  height  ZTP 
(cm) ;  default  =  48E5 


IV.  CALCULATION  OF  PRODUCTION  AND  LOSS  TERM 

The  coding  for  calculating  the  production  and  loss  rates 
is  set  up  in  JACOB.  For  example,  for  species  1^,  it  gives 


PR  ( I  x )  =  R(k1)SP(I2)SP(I3)  +  J(j1)SP(J4)  +  .  .  . 


1.0  ( 1 1 ) 


R  ( k .. )  SP  (Ir.)  +  J(j,)  +  •  .  . 


A-  14 


where  R(k^)  is  the  reaction  rate  for  reaction  k^, 

J(j^)  is  the  photolysis  rate  for  photolytic  reaction  3 ^ , 

SP(l2)  is  the  number  density  for  the  I2  species. 

With  all  values  corresponding  to  the  same  altitude.  The  K 

array  is  set  up  in  subroutine  COMR,  the  J  array  in  subroutine 

3f  i 

COMJ.  JACOB  also  give  the  Jacobian  —  necessary  for  solving 

dni 

the  photochemical  species. 

For  diffusive  species,  the  production  and  loss  rates 
are  stored  in  SSPC(81)  and  SSLC(81)  as  functions  of  altitude 
for  use  in  subroutine  DIFUSE. 


Calculation  of  Reaction  Rates  [COMR] 

Parameters  necessary  to  calculate  the  reaction  rate 
Rr(z)  are  stored  in  array  RR(I,K),  where 
I  =  reaction  parameter  index 

K  =  reaction  index  (see  BLOCK  DATA  BD2  for  list  of 
reactions) 

RR(1,K)  contains  the  parameter  defining  the  type  of 
reaction : 

a)  RR ( 1 , K)  =  1.0 
RK(z)  =  RR  ( 2  ,  K) 

b)  RR ( 1 , K)  =  2.0 

RR ( 3 , K) 

Rr(z)  =  RR(2,K)e  T(z)  ,  T(z)  -  tempera- 

cure  at 
altitude  z 


c)  RR  (1 ,  K) 
Rk(z)  = 


RR ( 3 , K) 

RR  ( 2  ,  K)  e~  T  (z  ) 
RR  (4 ,  K)  +  N 


t 


N  =  number  density 
of  air 


d)  RR(1,K) 


<1 


R  (z)  -  k°  Cl'UlMU)  o  6UiOy10^'^fT~('zf^'l2]  1 
RK(Z)  "  ,  .  k„  (T(z)M(z)  0,6  10  k°°  (TU) 

1  +  kM  (T  ( z )  ) 

where  k0(T(z)  =  RR (1 , K)  *  (y^-)  _KH  ( 2 '  k) 


kJCUzJU  RR(3,k)  *(|^-)"RR(4'k) 
The  parameter  RR(I,K)  can  be  changed  from  namelists. 


Calculation  of  Photolysis  Rate  [COMJ] 

The  photolysis  rate  J  is  defined  by 


J  (L)  =  1  QS(I)  xFLX(l)  xTAU(L,I) 

over  72  wavelengths 

t  h 

where  QS(I)  is  the  absorption  cross-section  at  the  I  wave¬ 
length 

FLX(l)  is  the  unattenuated  integrated  solar  flux  at 
t  H 

the  ]  wavelength  interval 

TAU  ( L ,  I )  is  the  transmission  function  at  the  L*"*1 

till 

altitude  and  IL  wavelength. 

The  transmission  function  is  calculated  in  the  subroutine 
TAU  where 

_ T 

TAU  =  e  cosO(t) 

with  t  the  optical  thickness  which  is  a  function  of  03  and 
ibution ;  0(t)  is  the  zenith  angle,  which  is  a  func¬ 
tion  of  latitude,  declination  and  local  time.  For  the  diurnal- 
average  option,  the  transmission  corresponds  to  the  averaged 
transmission  over  a  24  hr  day 
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TAU 


1  r  t/cosO(t) dt 
IThZJ  e 

over 
daylight 
hour 

V.  TREATMENT  OF  PHOTOCHEMICAL  SPECIES 

In  calculating  the  new  species  number  densities,  the 
user  has  the  option  of  either  solving  one  large  system, 
generated  by  considering  all  of  the  photochemical  species, 
simultaneously,  or  letting  the  program  group  the  photochem¬ 
ical  species  into  families  which  are  then  processed  sequen¬ 
tially,  as  sub-systems.  The  justification  for  this  approach 
comes  from  the  fact  that  the  coupling  among  the  species  with¬ 
in  each  family  is  much  stronger  than  the  coupling  to  species 
outside  the  family.  Presently,  there  are  four  photochemical 
families:  oxygen  (0<XD),  0) ,  hydrogen  (H,  OH,  H02,  h2o2, 
CH20)  ,  chlorine  (Cl,  CIO,  HC1 ,  H0C1 ,  ClNO-j)  and  nitrogen 
(NO,  N02,  HN03,  NO 3 ,  N205,  H02N02> . 

The  Diurnal- Averaged  Option 

In  the  diurnal-averaged  option,  one  has  the  set  of  equa¬ 
tions 


{P  (n)-L.  (n)xn.  =  0)  .  .  .  . 

i  i  =  l,m  (species), 


(A- 1 ) 


where  n  =  (n^,n2».  .  .nm)  are  the  species  densities, 
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Depending  on  whether  or  not  nil  of  the  photochemical  species 
are  being  solved  together  as  a  group  or  are  being  processed 
by  families,  the  number  densities  n  represent,  either  the 
entire  group  of  photochemical  species  or  just  one  particular 
photochemical  family. 

The  system  of  equation  (A-l)  is  solved  by  Newton's 
method  as  follows: 

Defining  the  set  of  functions 


(fi(n) 


P . 
1 


(n)  -Li  (ft)  x  n^} 


i  =  l,m 


(A- 2) 


and  expanding  the  f ^  in  Taylor  series  to  first  order  in  An,, 
about  n„  where  n0  is  either  a  guess  solution  or  a  solution 
from  the  previous  interation,  we  have 

m  i>f  .  (n0)rr  _ 

(£i(nl  -  V  0|4nl'  h-l.m 


where  An  .  =  n  .-n?  . 

1  3  3 

Since  f ^ (ft)  =  0,  for  i=l,  m,  the  following  system  of 
simultaneous  linear  equations  for  An_^  ,  the  first  order 
correction  to  n0 ,  is  obtained 


M  df i (n0) 


a: 

3  =  1 


9n  . 
3 


An. 


* 
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The  equations  are  generated  by  the  program  as  follows.  Using 
an  inital  guess  for  the  solution  n0 (supplied  by  either  BLOCK 
DATA  values  or  restart  values,  previously  calculated  and 
stored  on  disk),  subroutine  NEWTON  calculates  the  f^(n0)  and 


A- IB 


9fi <n0) 

subroutine  JACOB  calculates  the  Jacobian  terms  — ..  - 

J  n  • 

3 

For  the  chlorine  and  nitrogen  family,  the  system  of 
equations  (A-4)  cannot  be  solved  as  is  because  they  turn  out 
to  be  ill-conditioned.  One  can  define  the  species  concentra¬ 
tions 


C1X  =  Cl  +  CIO  +  HC1  +  H0C1  +  CINO3 

NOX  =  NO  +  N02  +  HN03  +  N03  +  N205  +  N205  +  H02N02 

+  cino3  . 

If  treated  as  species  on  their  own,  the  production  and  loss 
terms  for  C1X  (NOX)  are  just  the  sum  of  the  production  and 
loss  terms  of  each  of  the  species  on  the  right  hand  side. 

It  turns  out  that  this  net  production  and  loss  term  is  small 
compared  to  the  individual  contributions.  Thus,  the  terms 
on  the  right  hand  side  sum  approximately  to  zero  and  the 
system  of  equation  is  ill-conditioned  in  that  it  is  degener¬ 
ate.  To  remedy  this,  in  processing  the  chlorine  and  nitrogen 
family,  the  species  have  the  maximum  number  density  is  found. 

The  corresponding  equation  in  the  system  is  replaced  by  an¬ 
other  equation  formulated  from  number  density  closure  con¬ 
siderations.  In  the  case  of  chlorine,  we  require 

AnCl  +  AnC10  +  AnHCl  +  AnHOCl  *  AnClN03 

Clx  -  (Cl  +  CIO  +  HC1  +  H0C1  +  C1N03) 

In  the  case  of  nitrogen, 

SnN0  +  4nN02  +  4nHN03  *  4nN03  +  2  *4nN205  +  AnH0.,NO2  *  AnClN03 
=  N0x  -  (NO  +  N02  +  HN03  +  N03  +  2 x N205  +  H02N02  +  C1N03) 

is  the  replacement  equation.  The  two  equations  demand  that 
the  sum  of  the  odd  chlorine  and  nitrogen  species  is  conserved 
and  also  provide  some  feedback  from  diffusion  processes. 
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The  system  of  linear  equations  is  solved  in  FUNCTION 
LINEQN,  using  Gaussian  elimination  and  the  changes  in  den¬ 
sity  An.  are  returned.  The  new  densities  n  -  n0  +  An  are 
used  to  recalculate  the  production  and  loss  terms,  P^  and 

L  n  .  The  value  of  f.  =  P.-L.n.  is  then  found  and  the  ratios 
11  1111 

f./P.  are  tested  for  convergence  at  that  altitude.  If  | 
ii  ‘  i 

EP1,  convergence  has  been  achieved,  the  program  goes  on  to 
the  next  altitude.  Otherwise,  it  goes  back  and  recalculates 
all  of  the  families  again. 


Diurnal  Option 


In  calculating  the  photocnemical  species  under  the 
diurnal  option,  one  seeks  a  solution  to 


Dni (t) 


p.(n(tn  -  la  <H(t)  )xni(t)  )imltJn 

(species) 


(A- 5) 


with  a  periodic  boundary  condition  requiring  that  n^(t,z)  = 
n^ (t  +  24  hours,  z)  for  all  times  t. 

Equation  (A-5)  is  a  system  of  coupled  non-linear  differ¬ 
ential  equations  in  time.  One  seeks  the  solution  in  two 
steps.  First,  a  finite  difference  method  has  to  be  derived 


to  propagate  n(t)  in  time.  The  scheme  is  then  used  to  prop¬ 
agate  n(t)  in  time  until  the  periodic  condition  is  satisfied 
by  all  the  species.  A  periodic  solution  for  n(t)  exists  since 
the  driving  term  P^-L^n^(t)  has  the  24  hr  periodicity. 

The  finite  difference  scheme  is  first  discussed.  Using 
the  notation  that  n  represents  the  number  densities  at  the  k 
time  step,  we  can  rewrite  (A-5)  in  finite  difference  form 
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nk+1  -  nk 
t  Ti 


(At) 


k+1 


„  ,-k+l,  T  ,->k+l.  v  k+1 , 

P.(n  )  -  L.(n  )  x  n.  }i=lf 


m 


or 


(n1"1  -  At^lP^n"1)  -  L^n^*)  x  n~~l 


k+1 


k+1, 


k+1. 


(A- 6) 


-  n .  =  0}.  .  m 

i  i=l ,m 


k+1  k  A k+1 

where  (At)k+1  =  t  -  t  .  In  the  above  equation,  n  is 

-a  k 

the  unknown  and  it  is  assumed  that  n  has  been  computed  in 
the  previous  time  step.  We  define 


,  .j.k+1.  _  k+1  ..  ,.»k+l. 

(gi(n  )  =  n±  Atw,  IP,  (n  )  - 


'k+11  i 


T  ,-»k+l4  k+1,  k, 

L.  (n  )  x  n.  J  -  n.  }  .  ,  m 
i  i  i  i=l, m 


k+1 


(A-  7 ) 


The  required  solution  n  is  then  given  by 


(gi(*k+1>>i=1,In  -  o 


(A- 8) 


-*  k+ 1 

If  nQ  is  the  guess  solution  or  the  solution  from  the  pre¬ 
vious  iteration,  one  can  expand  g^  in  Taylor  series  to  first 
order  in  Ank+^  around  nk+1  to  obtain 


r  _  ,ik+l4  _  ,±K+i,  ,  „  'X  “  '  „_K+i, 

ig.  (n  )  2;  g.  (ne  )  +  E  - v-n —  An.  ).  , 

i  i  j=1  J  i"1' 


,*k+l.  ,  ”  agi(n°+1)  .k+1 


m 


(A-  9) 


k+1 

from  (A-8)  ,  (A-9)  gives  a  set  of  linear  equations  for  An..  , 
namely,  the  first  order  correction  term  to  nk+i, 
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(A- 10 ) 


, ;  3yftktli  * 

3  =  1  3  J 


,->-k+l.  . 

_gi (n°  i=l  ,m 


Note  that 


agi(hk+1) 

9nk+1 

3 


3f i (nk+1) 

’’an*'*"" 


1- At 


Of. (nk+1 ) 

-SF- 

3 


if  i  /  j 


if  i  =  j  , 


„  ,-»k+l . 

Og . (n  ) 

Thus,  subroutine  JACOB  can  be  used  to  generate  — 1^+j- - 

9n  j 

with  minor  modification.  The  system  of  linear  equations 
( A- 1 0 )  is  solved  using  subroutine  LINEQN  to  return  the  den- 
sity  changes  An^  .  The  new  densities  n  =  nD  +  An 
are  used  to  recalculate  the  production  and  loss  terms  P^  and 
Lini‘  The  gi  are  then  f°und  and,  after  being  normalized, 
are  tested  for  convergence  against  EP1.  If  convergence  has 
been  achieved,  the  program  goes  on  to  the  next  time  step. 

If  not,  it  goes  back  and  recalculates  all  of  the  fast  families 
again . 

The  equations  are  then  propagated  forward  in  time  until 
the  24  hour  periodic  boundary  conditions  are  satisfied.  Sub¬ 
routine  SETUP  is  then  used  to  calculate  diurnally-averaged 
fast  species  densities  and  photolysis  rates,  as  well  as  cer¬ 
tain  terms  which  contribute  to  the  production  and  loss  of 
selected  slow  species.  These  values  are  saved  in  prepara¬ 
tion  for  entering  the  slow  chemistry  phase  of  the  calculation. 
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VI. 


DIFFUSIVE  SPECIES 


I 


\ 


The  diffusive  species  are  comprised  of  N2Of  ,  n2, 
ch-jCI,  cci4,  cfci3,  cf2ci2/  ch3cci3,  cf2o,  CO,  XX,  nox,  cix 
and  03 .  The  slow  species  concentrations  are  obtained  by 
solving  the  diffusion  equation  one  species  at  a  time.  The 
diffusion  equation  can  be  written  in  the  form 


^(KN||)  =  P  (z)  -^(z)f(z) 


( A- 11 ) 


where  z  is  geometrical  altitude,  K  is  eddy  diffusion  coeffi¬ 
cient,  N  is  air  density,  P  is  the  product  rate  and  (z)  is 
the  loss  frequency  given  by  LxN. 

The  differential  operator  is  converted  to  a  finite  dif¬ 
ference  operator  using  central  differencing. 

Using  the  notation  that  f^  stands  for  f(z^),  we  have 


s,KNaf' 


(KN)  .  -  _ 


Az 


Az 


1KN)  j-* 


or 


i 


- 


(A- 12) 


(KN)i-tsfi-1  +  ('(KN)  -  (KN)  ,+}s}f.  +  (KN)  .^f.+i 

(Az)2 


Thus,  ( A- 1 1 )  reduces  to 
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-  — ^  - 


and  F .  =  -P .  .  Note  that  (A-13)  holds  only  for  j=2 
D  D 

through  n-1  where  n  is  the  top  altitude. 

To  find  the  production  and  loss  terms  P^  and the 
program  calls  FUNCTION  FUNPL.  Calculation  of  the  eddy  dif¬ 
fusion  coefficients  K ^  is  done  in  subrouting  FUNK  and  the 

A.,  B . ,  C.,  F.  are  calculated  in  subroutine  ABCF . 
Dili 

The  above  set  of  equations  will  generate  a  system  of 
n-2  linear  equations  in  the  n  unknonws  f...  In  order  to  pro¬ 
duce  two  additonal  equations,  bottom  and  top  boundary  condi¬ 
tions  arc  used.  Either  mixing  ratio  or  flux  boundary  condi¬ 
tions  may  be  used.  When  mixing  ratio  boundary  conditions 
are  chosen, 


f ^  =  given  mixing  ratio  of  the  species  at  altitude  z, 
and/or 

f n  =  given  mixing  ratio  of  the  species  at  altitude  zn 
are  the  equations  provided.  If  flux  boundary  conditions  are 

j 

used,  we  have  | 


1 


--2 (KN)  - 

f  4- 

■2 (KN) 

.  ( Az ) 2  . 

.(Az)2  . 

-44> 


1 


0.5 


(KN) 1i 
<KN)  " 


/  (2Az) 
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" (KN)  .  +  (KN)  1 
n-1  n 

•KN»n-l  +  (KNI„" 

2 

V.  ^  -L 

L  2  J 

AZ  n-1 

Az 

(A- 16 ) 


=  4> 

n 


where  <t> ^  and  <J>n  are  the  specified  flux  at  the  bottom  and  top 
respectively . 

Hence,  including  the  two  equations  provided  by  the  bot¬ 
tom  and  top  boundary  conditions,  we  have  a  system  of  n  linear 
equations  in  the  n  unknowns  f^.  Furthermore,  the  system  is 
tridiagonal  and  is  solved  in  FUNCTION  TRIDIA,  which  returns 
the  new  mixing  ratios  at  altitudes  z^_^  n* 

As  each  slow  chemistry  species  is  calculated,  the  new 
mixing  ratios  f ^  are  converted  back  to  number  densities  and 
compared  with  those  from  last  iteration.  If  the  relative 
changes  are  all  less  than  a  pre-specif ied  value,  at  all  alti¬ 
tudes,  then  that  species  is  determined  to  have  an  acceptable 
density  profile.  If,  on  a  particular  interation,  all  slow 
species  have  acceptable  new  density  profiles,  then  the  entire 
set  of  species,  both  fast  and  slow,  is  considered  to  have  con¬ 
verged.  Otherwise,  the  program  recalculates  the  O^  column 
densities  and  then  the  photolysis  rates  in  beginning  the  next 
iteration  and  recalculation  of  the  fast  chemistry. 
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VII  OUTPUT 


In  the  diurnal-averayed  option,  the  level  of  output  is 
conrolled  by  the  parameter  LEVP 

LEVP  =  1  PRINT  altitude  array  for  densities  and  mixiny 
ratio  for  diffusive  species 

LEVP  =  2  PRINT  altitude  array  for  photochemical  species 

LEVP  =  3  PRINT  altitude  array  for  fixed  species  (M, 

H2O,  O2) /  temperature,  eddy  diffusion  coeffi¬ 
cients  and  O2  column  density 
LEVP  =  4  PRINT  reaction  rate  tables  and  boundary  con¬ 
ditions 

LEVP  =  5  PRINT  altitude  array  for  photolysis  rates 
and  column  density 

LEVP  =  6  PRINT  selected  reaction  rates  and  calculated 
fluxes 

LEVP  =  7  PRINT  wavelength,  solar  flux  and  absorption 
cross-section . 

The  program  is  set  up  so  that  all  options  lower  than  '.-.ho 
specified  number  are  included  in  the  print  out.  In  addition, 
the  altitude  arrays  are  controlled  by  che  parameters 
ZBP(cm):  bottom  level  to  be  printed  out 

ZTP(cm) :  top  level  to  be  printed  out 

UZ(em):  increments  between  orint  out  levels. 

A  listing  of  a  sample  output  is  included  in  this  Appendix. 
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VIII . 


RUNNING  INSTRUCTIONS 


The  program  is  currently  being  run  on  the  AFGL  CDC  6600. 
A  typical  deck  setup  together  with  a  typical  namelist  input 
is  given  in  Table  a-4. 

The  typical  running  time  of  the  diurnal  average  version 
is  about  60-80  seconds  on  the  CDC  6600.  For  the  diurnal 
version,  the  time  required  is  approximately  2,000  sec.  for 
a  complete  run.  If  one  only  executes  the  part  to  solve  the 
photochemical  species  only,  the  time  required  is  reduced  to 
about  600  seconds. 
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Table  A- 4 


PROGA,  CM120000,  T100-  1234 

FTN ,  I  =  INPUT,  OPT  =  2,  SL,  R  =  3. 

COPYCF ,  INPUT,  TAPE 5 . 

REWIND,  TAPE 5 . 

COPYCF,  TAPE 5 ,  OUTPUT. 

REWIND,  TAPE5 . 

LDSET  (PRESET  =  ZERO). 

LGO . 

7/8/9 

FORTRAN  DECK  (Fortran  IV) 


7/8/9 
2  L 1 3T 

FSF  (1  ')  =,'  .  , 

ICF  (  5 ,1  )=  1  ,  IGF ( 3, Z)  -  >1  , 

ICF  ('♦ ,  1 )  =  2  »  ICF  ( N ,  2  )  =5  Q„  TCP  ( ,  v  )  - 1  2  , 
RR  (  1  ,'?)=<?.  6-3-3  C  , 

R  R  ( 1 ,  - 1 )  ~  c  .  1 E  -  3 1  , 

F JC  3=  2. 

R  C  u  •  t 

k=i 

mi  v=a, 

pctz= : , rrc=5Q , 

EP1=2*l'  E-4»EP  2"=  £•.,£“  if 
ZE-C.  1,  /T  =7  0.  tit  5, 

2EP  =  J  ,J,Z  f  P  =  7  0 »  L  L  3  »  C  Z  -  1  •  £  » 

IE  VP=  TRAN=  .  TPl  F  .  , 

E 1 1  =1  ,i.-  it  H2  =  .  3  , 

E C B  <19  5  =  3 .t-7, 

CCS (3 J) -1 . 5E-F, 

E  C  P  ( ?.  1 1  --  9  .  C  -  7  , 
eCE  (22)  =1  .2F-9, 

EC  n  <  I  .')  .  30  1C  , 

ecn  <0)  =•!  .  30  9, 

FCR  < 1  r>)  r  j  ,  v-  +  f  , 

ECS  (  ;1  -•>  =!.•:/» 

E T \ 2  7)  =1  .  t  ♦  3, 
eco  (.?  =  3 .o 7, 

EC”  (  7  )5  -  t  .  O  3 , 

6  C  f<  (  3 1  J  -  1  .  t  -  1  J  j 

ec« ( 1 1) -i ,c- ) , 

OCR  (  57)  -1  .  rt  C  -  i  E-  0  T  (  3  2)0.  j  f 
SSI  t  -»T’  .  , 
it  no 
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OPTION 


Figure  A-l  (cont.) 
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Calculate  Production  and  Loss 
Ratos  for  Diffusive  Species 

using  Exact  Diurnal  Averaging 

- 


Figure  A-l  (cont.) 


Solve  Diffusion 
Equation  for 
Diffusive  Species 
One  at  a  Time 


Update  all 
Number  Densities 


Odd  Ntrogcn 
Densities  so 
that  they  are 
Consistent  with 
Newly  Calculated 
N0X  and  Clx 


Figure  A- 1  (cont.) 


APPENDIX  B 


DESCRIPTION  OF  2-D  MODEL 


I.  General  Background 

II.  Chemical  Scheme 

III.  Treatment  of  Photochemical  Species 

IV.  Treatment  of  Diffusive  Species 
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1.  GENERAL  BACKGROUND 


The  atmosphere  is  divided  into  grids.  The  coordinates 
used  are  the  latitude  <p  and  log-pressure  coordinates 


£  =  In  ff. 

P 

where  p0  is  the  surface  pressure.  In  this  coordinate,  the 
vertical  velocity  w  =  ^  is  related  to  w  =  ^  via 


3  _  3 

^  "  3€ 

Using  the  coordination  transformation ,  we  can  rewrite  (7.2-10) 
and  (7.2-14)  as  follows 


,  -L  V  .  c  \  -L  U 

3t  a  cos(i  3A  u  a  cosiJ>  3$ 

r  .3  .  r  -  r,  , 

e  e  w)  =  Q/p 

1  3u  1  3  .  .. 

acos(i  3A  a  cosi)  34>  v  c°s<J>) 


(  f  V  COS<)>)  + 


+  •—  (e~f'w) 


=  0 


(B-l) 

(B  —  2 ) 


Introducing  the  variables 


B-2 


V 


(B-3) 


=  v  e  ^  cos4> 
W  =  u>  e  ^  cost p 


one  can  recast  (B-l)  and  (b-2)  in  the  form 


3f 


(fu)  + 


es  3 


3t  a  cos4>  3  A '  a  cos<J>  3<f> 


(£  V)  + 


e.g  3  =  Q 

0os<fr  dr.  '  p 


3u  *  -Aiv  + 


a  cos<t>  3  A  acosiji  34>  cos<{>  3£ 


=  0 


Upon  zonal  averaging,  one  obtains 


(B-4 ) 


(B  -  5 ) 


3  f 

e^ 

-1 

3 

at 

cos<f> 

a 

3  (j> 

13V 

3W 

0 

a  3$ 

9C 

(vf)-  A(wf) 


Fr  +  —  e  ^  cosi| 
f  P 
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(B  -  7 ) 


Equation  (A-7)  implies  the  existence  of  a  stream  function  4*, 
where 


V 


3  T 


w 


1  34' 
a  3$ 


(B  -  8 ) 


Finally,  the  eddy  fluxes  are  parameterized  by  eddy  diffusion 
tensor,  giving 


B-J 


i 


3 

34 


(B-9) 


i r  -  4 

a  3  <J>  e 


cos*l%*S  +  Noe' 


-4  ,  3f  „  3f  . 

e  cos>i>(K.  .  .  +  K.r  ) 

4<J>  3<J>  44  3t, 


The  globe  is  covered  by  19  latitude  belts  each  '9.5° 


wide  (or  '^1,000  km  across).  The  vertical  direction  is  di¬ 


vided  into  layers  with  equal  spacing  in  4  with  A4  -  0.5  (or 


'^3  km  high)  .  The  model  is  set  up  to  accommodate  a  maximum 


ol"  2  9  layers. 


Fig.  B-i  shows  a  schematic  diagram  of  the  grid  together 


with  typical  points  within  a  grid  for  which  variables  are 


defined.  Variables  are  defined  at  points  A  (at  center  of 


boxes)  with  indices  J,  K  with  the  corresponding  coordinates 


given  by 


.  /  r  1  \  TT  .  IT  11 

4>  (J  1)  X  19  +  38  2 


J  =  1,19 


r  ~ 


+  0.25 


K  =  1,29 


SP ( I , J, K) :  diffusive  species  mixing  ratio,  I  is  species 


index 


T  ( J,  K) 


temperature 


M(J,K) : 


air  number  density 


DADT ( J , K ) :  local  time  rate  of  change  of  mixing  ratio 


Variables  are  defined  at  points  C  and  G  with 


=  (J-i)  +  £  -  J 


J  «  1,19 


c  = 


K  =  1,30 


I 


K5C(JfK)  : 
W(J,K) : 

VF  ( J ,  K)  : 


component  of  eddy  diffusion  tensor 
zonal  mean  vertical  wind  speed  W 
vertical  eddy  fluxes  (Kf ~  ijy-) 

Variables  are  defined  at  points  B,  D,  F,  H  with 

4>  =  (J-l)  Y9  ~  ^  J  =  1,20 

f  “  K  =  1,30 

T(J,K):  stream  functions  that  general  V,  W 

K^_(J,K):  diagonal  element  of  eddy  diffusion  tensor. 

Variables  are  defined  at  points  E,  I  with 

«t»  -  (J-l)  x  ^  -  ~  J  =  1,20 

f  +  0.25  K  =  1,29  . 

K^(J,K):  triagonal  component  of  eddy  diffusion  co¬ 

efficient 

V ( J , K) :  zonal  mean  horizontal  velocity  V 

HF  ( J ,  K)  :  horizontal  eddy  fluxes  (K,  .  ^  +  K..  -jr  )  . 

4><t>  <K 

The  variables  are  defined  in  this  way  so  that  upon  applying 
the  finite  difference  scheme  to  obtain  the  appropriate  spatial 
derivatives,  all  quantities  that  enter  into  (R-6)  are  deiined 
at  the  centers  of  the  boxes. 

The  overall  proceeding  for  solving  the  atmospheric  species 
is  outlined  in  Fig.  B-2.  At  present,  all  the  dynamic  variables 
are  read  from  a  local  file.  The  temperature,  stream  function 
and  the  photochemical  species  concentrations  are  updated  every 
10  days,  while  the  eddy  diffusion  tensor  is  updated  once  every 
30  days.  Some  of  the  components  outlined  in  Fig.  b-2  will  be 


discussed  in  the  following  sections. 


II.  CHEMICAL  SCHEME 


The  chemical  scheme  in  the  2-D  model  is  similar  to  that 
set  up  in  the  one-dimensional  model.  The  species  are  divided 
into  the  photochemical  species,  diffusive  species  and  the 
fixed  species. 

The  photochemical  species  include: 

Oxygen  family:  0(^D),  0(^P) 

Hydrogen  family:  H,  OH,  H02,  H2°2 

Chlorine  family:  Cl,  CIO,  HC1,  HOC1 ,  C1NO 3 


Nitrogen  family:  NO,  N02, 


HNO- 


The  diffusive  species  are: 


CH3C1,  CH3CC13.  CC14,  F-ll,  F- 12 ,  C1X,  N20,  NOX 
and  03 

The  fixed  species  are: 

ch4,  CO,  h2o. 

The  chemical  reactions  and  photo  lytic  processes  adopted 
correspond  to  the  appropriate  subset  of  the  1-D  model  chem¬ 
istry  as  listed  in  Appendix  A.  The  boundary  conditions  are 
also  similar  except  for  F-ll  and  F-12.  A  time-dependent  flux 
boundary  condition  is  used  for  both  species.  The  release  of 
F-ll  and  F-12  is  restricted  to  between  25°  N  and  60°  N.  The 
release  rate  follows  the  data  compiled  by  the  CMA  Fluorocarbon 
Panel  (1980).  The  values  are  given  in  Table  B-l  for  reference. 
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TABLE  B-l 


FLUOROCARBON  RELEASE  RATES 

FLUCROr-'P.iiC.'l  r5;rnn,~ri0..,  AK0  RELEASE 
iULiION  KIEOOAAMa 

REFORTItw  CO.IPAHIES  ONLY 


IC-11 


AliUUAL  TOTAL 


YEAR 

Pa  CD 

UEL* 

p;;oo 

HCL 

UNUCL’ 

19f>0 

6.6 

5.6* 

10.6 

14.4 

4.1* 

1951 

9.  1 

7.5’ 

27.6 

21.9 

5.  7' 

1932 

13.6 

10. 6‘ 

41.2 

32.7 

0.5' 

1553 

17.3 

14.7' 

53.5 

*7.4 

11.1' 

USA 

20.9 

13.3' 

79.'* 

65.7 

13.7' 

1955 

C6 . 3 

22.6' 

105.6 

08.3 

17.4' 

1956 

32.5 

28.2* 

133.1 

116.5- 

21.6' 

1  55  7 

33.9 

31.6' 

172.0 

143. 1 

2  3.9* 

1953 

29. 5 

29.7' 

201.6 

177.8 

23.6' 

1959 

35.6 

30.3' 

237.1 

203.1 

29.0' 

lV'O 

V).  7 

39.7* 

CO 6. 9 

247.9 

39.0* 

1961 

oO .  5 

51.2! 

34  7.3 

259.0 

48.3' 

1962 

73.1 

«>*»  .  1  ' 

625.4 

363.1 

62.3' 

1961 

93.3 

73 . 5  ' 

516.7 

441.7 

77.0' 

1969 

111.1 

93.2' 

629.3 

534.9 

94.9' 

196  5 

ICC.  a 

106.3' 

7 : » C  .  6 

641.2 

111.4* 

1966 

191.0 

119.0' 

893.7 

760.2 

133.4' 

156  7 

159. S 

135.1' 

1053.4 

695.3 

153.1' 

1963 

103. 1 

153.9' 

1236.5 

1049.2 

’87.4' 

1969 

217.3 

173.4' 

1453.8 

1227.6 

226.2' 

1970 

233.1 

202 . 8 ' 

1691.9 

1430.4 

261.o' 

1971 

263.2 

2CC . 7* 

1955.1 

1653.1 

302.0' 

1972 

306.9 

250.0' 

2Co2 . 0 

1903.9 

358.0' 

1975 

35  9. 1 

236.8' 

2611.1 

2190.7 

420.5' 

1979 

369.7 

315.4' 

2530.8 

2506.1 

474.7' 

31** .  1 

52 94 . 0 

?Z\  l  .6 

4  S  3  .  P  ' 

1976 

339.3 

295. S' 

2o34 . 7 

3107  .4 

527.3' 

19 ;/ 

320.5 

297.2' 

3955. ! 

3404 . 6 

550.6' 

19/6 

505.9 

231.0', 

4Co4 . 0 

io .  0 

578.4 ' 

1979 

269.5 

161.9' 

4553.5 

3947.5 

606.0’ 

AtirrjAL  ‘  TO  [  At 


YEAH 

PPOD 

PL  L  * 

P.iOU 

KLL 

UKWlL* 

1950 

34 . 0 

27.1' 

193.1 

129.6 

60.0' 

1951 

36.2 

30.2' 

234.5 

159.6 

7: .  9  ’ 

1952 

37.2 

31.5' 

2/1.7 

1«1  .  1 

60. 0’ 

1°53 

46.5 

35.5' 

316.2 

22o  .6 

91.6' 

1554 

49.1 

40. 3 ' 

3o  7 . 4 

2u6 . 9 

100.4' 

1 555 

57.6 

45.2' 

425.0 

312.1 

112. 9 ' 

1936 

63.7 

52.6' 

493.6 

36*.  7 

129.0* 

1  c  5  7 

74 . 2 

69.  .5* 

567.0 

4  24 . 6 

14  3.1' 

1930 

73.4 

62.6* 

641.2 

467.1 

154,2' 

1539 

37. 6 

6°.  fa' 

728.8 

'  5  ”-,  7 

172.1' 

lco0 

99.4 

1. 5 .  ' 

020  .  1 

1, . ') 

188. 3' 

1961 

103.5 

93.2' 

93o.  3 

733.2 

203.6' 

1562 

123.1 

'07.1' 

10c-4. 9 

840.3 

224.6' 

19o3 

1*6.4 

125.0' 

1211.3 

910.1 

2*6.1' 

1964 

170.1 

146 . 6 ’ 

1381 .4 

1112.7 

263.7' 

14*5 

160.1 

166.0* 

1 671.4 

i.:.3.4 

243.1 ' 

1966 

216.2 

184,3' 

1/07.6 

14c  2 . 7 

324.9' 

19f  7 

242.8 

208.3' 

2030.4 

16  70.9 

365.4' 

19c-S 

2  0  7 . 5 

233.0' 

2297.9 

150  4 .6 

343.1 ' 

1969 

297.3 

260.1’ 

25'>5.1 

ciG*  .6 

430.5' 

1570 

321.1 

234.2' 

2916.2 

24*3.9 

4  6  /  .  4  ' 

1971 

341 .6 

304.0-' 

3257.3 

2753.7 

804.1' 

1972 

379.9 

331.2* 

3637. 7 

3034 . 9 

552.6' 

1973 

423.3 

366 . 5 ' 

4081.0 

3451.4 

609.6' 

1974 

442.6 

396.5* 

4503.3 

33  .7.9 

655.9' 

1976 

331.0 

3 92.0' 

4304  .  Y 

42 10.0 

6 1 . 4 . 0  ' 

1976 

410.7 

3o3 . 9 ' 

S295.5 

4554 .6 

700.9' 

14  77 

3.52 . 0 

344 . 9 ' 

6473. 3 

4 4  A 9.0 

7  3."../)  • 

lWd 

372.1 

5io . 2  * 

o050.4 

6256.0 

/V*  .6' 

1979 

357.2 

305.0’ 

6  4  0  7 . 0 

55;1.6 

64a. O' 
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ill.  TREATMENT  OF  PHOTOCHEMICAL  SPECIES 

The  treatment  of  photochemical  species  is  exactly  anal¬ 
ogous  to  the  diurnal-averaged  version  of  the  1-D  model. 

Since  there  is  no  spatial  coupling  among  the  equations,  the 
photochemical  species  can  be  solved  locally  and  there  is  no 
difference  between  the  1-D  and  2-D  models.  A  detailed  de¬ 
scription  of  the  method  can  be  found  in  the  corresponding 
section  in  Appendix  A. 

IV.  TREATMENT  OF  DIFFUSIVE  SPECIES 


An  explicit  time  difference  scheme  is  used  to  propagate 
the  specie  volume  mixing  ratio  forward  in  time.  Given  the 
current  values  of  f,  v,  w,  K  and  T  at  time  tk,  the  program 


calculates  each  individual  term  on  the  right  hand  side  of 


(F-6)  to  obtain  (tj— )  ^ 


The  value  of  f  at  time  tk+1  is 


obtained  by  the  Adam-Bashford  Scheme, 


'tk+1 


i  5  (Uq 
•  'at* tk 


4t  -  ° - 5  >  tk-i At 


(B-10 ) 


Because  an  explicit  time  scheme  it.  used,  the  time  step 
must  be  sufficiently  small  to  ensure  stability.  With  the 
present  choice  of  grid  sizes,  a  time  step  of  4  hours  is  re¬ 
quired  . 


1,000  mb  i  vl^ 


90°  S 


90° 


N 


J,  latitude  (19  boxes) 


Figure  B-l 

Diagram  Showing  the  Grid  System  Employed 
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Figure  B-2 


